Principles of the polynomial Multiplication Algorithm and corresponding C code implementation --- Using FFT

Recently, I have carefully studied the quick computing of Polynomial multiplication in the introduction to algorithms, mainly using FFT, which is also implemented by myself. The summary is as follows.

1. polynomial Multiplication

Two polynomials are multiplied by polynomials, for example: 3 * x ^ 7 + 4 * x ^ 5 + 1 * x ^ 2 + 5 and 8 * x ^ 6 + 7 * x ^ 4 + 6 * x ^ 3 + 9 submultiplication, A Polynomial with a maximum number of values of 13 is obtained. In general, the general calculation method is: multiply each item in polynomial A and each item in polynomial B to obtain n polynomials, and then add each polynomial together, to obtain the final result, assuming that the maximum number of times a and B is n-1 and the length is N, the final calculation result requires O (N ^ 2) time complexity. However, using Fast Fourier Transform (FFT) can reduce the time complexity to O (nlog N ). This is because the time complexity of positive/inverse fast Fourier transformation for a complex number of sequences is O (nlog N ), the transformed sequence is a result of Polynomial multiplication for the original sequence (polynomial multiplication is equivalent to Convolution ). Therefore, FFT can reduce the time complexity of Polynomial multiplication. The specific explanations and proofs are described in the Introduction to algorithms or any other related algorithm book. In addition, it should be noted that some other operations can also be converted into polynomial multiplication, and then the FFT can be used to speed up the operation. Example: 1. the number multiplication operation is similar to the polynomial multiplication. The operation of a * B is to multiply the number on each digit of A by the number on each digit of B. FFT can greatly speed up operations, especially in the multiplication of large numbers. 2. For detailed data of different length paths from A to B, and detailed data of different paths from B to C, calculate the number of different length paths from A to C. We can regard the paths from A to B and from B to C as different times. For example, there are three paths from A to B with a length of 4 or 2 with a length of 5, B to C has one path with 2 or 4 lengths and 3, so the number of paths with different lengths from A to C is equal to (3 * x ^ 4 + 2 * x ^ 5) * (4 * x ^ 3 + 1 * x ^ 2) after the coefficients are converted to the polynomial multiplication problem, FFT can be used to speed up the computation.

2. FFT

Most people only need to use FFT, but this algorithm is quite basic. So I programmed and implemented it myself. The total code is only about 150 lines, which is not long. Of course, if the input sequence length is not an integer power of 2, I have no corresponding preprocessing, which is a lazy. In fact, you only need to take the logarithm of the length. For example, if the input length is 37, first expand 37*2 to 74 (which is required by FFT ), then, you can perform an integer operation on 74 log2 to get 27 = 128. Therefore, add 54 more zeros after 74.

In addition, it should be noted that the reverse function is required in FFT and IFFT. However, given that polynomial multiplication requires paired FFT and IFFT, during polynomial multiplication, reverse can be omitted (of course, this function takes a small amount of time ). The value of W should be calculated in advance. In this way, when calculating each item in FFT and IFFT, W does not need to be repeatedly calculated, which can save a lot of time.

The principles and explanations of FFT can be used to look up textbooks in any field, such as Wikipedia, information theory, digital signal processing, and random processes.

3. Code

The program mainly includes FFT, IFFT functions, and some complex operations.

3.1 Definitions of plural numbers and related operation Definitions

// Complex struct Complex {double real; double image ;}; complex A1 [max_size], A2 [max_size], result [max_size], W [max_size]; // calculate complex operator * (complex A, complex B) {complex R; R. real =. real * B. real-a.image * B. image; R. image =. real * B. image +. image * B. real; return r;} // calculate complex operator + (complex A, complex B) {complex R; R. real =. real + B. real; R. image =. image + B. image; return r ;}// complex subtraction complex operator-(Co Mplex A, complex B) {complex R; R. real =. real-b.real; R. image =. image-b.image; return r;} // complex division complex operator/(complex A, double B) {complex R; R. real =. real/B; R. image =. image/B; return r;} // complex operator ~ (Complex a) {complex R; R. Real = A. Real; R. Image = 0-a.image; return r ;}

3.2fft and IFFT functions and related functions

In fact, the principle of FFT is the same as that of IFFT, except that IFFT has an additional Division step and can combine the two into a function. Reverse is used to rearrange the element subscript of the input array. For example, if the input array is 8 characters in length, the elements of the following indexes become, the elements of the base object. Compute_w is used to calculate the value of W in FFT in advance.

// Rearrange method 2 with High Efficiency void reverse (int * ID, int size, int m) {for (INT I = 0; I <size; I ++) {for (Int J = 0; j <(m + 1)/2; j ++) {int V1 = (1 <(j) & I) <(m-2 * J-1); int v2 = (1 <(m-j-1) & I)> (m-2 * J-1 ); id [I] | = (V1 | V2 );}}};

// Rearrange method 1, which is less efficient by using the POW function void reverse (int * ID, int size, int m) {for (INT I = 0; I <size; I ++) {for (Int J = 0; j <m; j ++) {int exp = (I> J) & 1; id [I] + = exp * (INT) Pow (double) 2, (double) (m-j-1 ));}}};
// Calculate and store the required W value void compute_w (complex W [], int size) {for (INT I = 0; I <size/2; I ++) {W [I]. real = cos (2 * pI * I/size); W [I]. image = sin (2 * pI * I/size); W [I + size/2]. real = 0-W [I]. real; W [I + size/2]. image = 0-W [I]. image ;}}; // fast Fourier void FFT (complex in [], int size) {int * id = new int [size]; memset (ID, 0, sizeof (INT) * size); int M = Log (double) size)/log (double) 2); reverse (ID, size, M ); // re-arrange the input to conform to the output complex * resort = new complex [siz E]; memset (resort, 0, sizeof (complex) * size); int I, J, K, s; for (I = 0; I <size; I ++) resort [I] = in [ID [I]; for (I = 1; I <= m; I ++) {S = (INT) Pow (double) 2, (double) I); For (j = 0; j <size/s; j ++) {for (k = J * s; k <j * s + S/2; k ++) {complex k1 = resort [k] + W [size/S * (k-J * s)] * resort [K + S/2]; resort [K + S/2] = resort [k]-W [size/S * (k-J * s)] * resort [K + S/2]; resort [k] = k1 ;}}for (I = 0; I <size; I ++) in [I] = resort [I]; Delete [] ID; Delete [] resort ;};// Fast inverse Fourier void IFFT (complex in [], int size) {int * id = new int [size]; memset (ID, 0, sizeof (INT) * size ); int M = Log (double) size)/log (double) 2); reverse (ID, size, m); // rearrange input, complex * resort = new complex [size]; memset (resort, 0, sizeof (complex) * size); int I, J, K, S; for (I = 0; I <size; I ++) resort [I] = in [ID [I]; for (I = 1; I <= m; I ++) {S = (INT) Pow (double) 2, (double) I); For (j = 0; j <size/s; j ++) {for (k = J * s; k <j * s + S/2; k ++) {C Omplex k1 = (resort [k] + (~ W [size/S * (k-J * s)]) * resort [K + S/2]); resort [K + S/2] = (resort [k]-(~ W [size/S * (k-J * s)]) * resort [K + S/2]); resort [k] = k1 ;}}} for (I = 0; I <size; I ++) in [I] = resort [I]/size; Delete [] ID; Delete [] resort ;};

3.3 Main Function

Input the coefficients of two polynomials (the length must be an integer power of 2) and output the result of multiplying two polynomials.

Int main () {// input two Polynomial Series int size, size1, size2, I; memset (A1, 0, sizeof (A1); memset (A2, 0, sizeof (A2); memset (W, 0, sizeof (w); memset (result, 0, sizeof (result); scanf ("% d ", & size1, & size2); for (I = 0; I <size1; I ++) scanf ("% lf", & A1 [I]. real); for (I = 0; I <size2; I ++) scanf ("% lf", & A2 [I]. real); size = size1> size2? Size1 * 2: size2 * 2; compute_w (W, size); FFT (A1, size); FFT (A2, size); for (I = 0; I <size; I ++) result [I] = A1 [I] * A2 [I]; IFFT (result, size); for (I = 0; I <size1 + size2-1; I ++) printf ("%. 2lf ", result [I]. real); printf ("\ n"); Return 0 ;}

The complete code is as follows:

CPP File Download

Principles of the polynomial Multiplication Algorithm and corresponding C code implementation --- Using FFT