polynomial multiplication
given a $n$ degree polynomial $A(x)$ and $m$ degree polynomial $B(n)$ , find out $F(x)$ with $G(x)$ of convolution.

Input format

Two integers in the first line $n,mn,m$ 。
next line $n+1$ numbers, expressed from low to high $F(x)$ coefficient.
next line $m+1$ numbers, expressed from low to high $G(x)$ coefficient.

Output format

one line $n+m+1$ numbers, expressed from low to high $F(x) dot G(x)$ coefficient.

Input and output samples

Input sample

1 2
1 2
1 2 1

Output sample

1 4 5 2

Problem analysis

Suppose two polynomials $A(x)=a_0+a_1x+a_2x^2+dots +a_{n-1}x_{n-1}$ and $B(x)=b_0+b_1x+b_2x^2+dots +b_{n-1}x^{n-1}$ , the two polynomials can be written as:
$egin{equation} A(x)=um_{i=0}{n-1}a_ixi B(x)=um_{i=0}{n-1}a_ixi nd{equation}$
The traditional method is to use the coefficients of two polynomials to perform a convolution operation to obtain a $2n-2$ degree polynomial $C(x)$ ：
$egin{equation} C(x)=um_{i=0}^{2n-2}c_ix_i nd{equation}$
The time complexity of this convolution operation is $O(n^2)$ , obviously it is unbearable when the data range is large, and usingFast Fourier TransformThe time complexity can be reduced to $O(nlogn)$ 。

Introduction to FFT

Fast Fourier transform (fast Fourier transform) is the collective name for efficient and fast calculation methods that use computers to calculate discrete Fourier transforms (DFT), referred to as FFT. Fast Fourier transform was proposed by J.W. Cooley and T.W. Tukey in 1965. Using this algorithm can greatly reduce the number of multiplications required by the computer to calculate the discrete Fourier transform. In particular, the more sampling points N are transformed, the more significant the calculation savings of the FFT algorithm will be.
FFT (Fast Fourier Transformation) is a fast algorithm for Discrete Fourier Transform (DFT). That is the fast Fourier transform. It is obtained by improving the discrete Fourier transform algorithm based on the odd, even, imaginary, real and other characteristics of the discrete Fourier transform.
——Baidu Encyclopedia

The problem to be solved is the multiplication of polynomials, and the representation method of a polynomial is not unique. Traditionally, the representation of polynomials uses the coefficient representation method, that is, for a polynomial $A(x)=a_0+a_1x+a_2x^2+dots +a_{n-1}x_{n-1}$ can be determined by a coefficient vector $(a_0,a_1,a_2,dots,a_{n-1})$ The only representation.

In addition to coefficient representation, polynomials can also be represented by point values. For polynomials $A(x)$ , selected $n$ a $x$ value $(x_0,x_1,x_2,dots,x_{n-1})$ Bringing in polynomials for calculation, we get $n$ pip value) this $n$ Point values can also uniquely represent polynomials $A(x)$ 。

So when we use the same vector $x$ After obtaining the point value representations of two polynomials, multiply them by the corresponding point values to get $(A(x_0)B(x_0),A(x_1)B(x_1),A(x_2)B(x_2),dots,A(x_{n-1})B(x_{n-1}))$ That is, the point value representation of the product polynomial of two polynomials. The time complexity of this process is $O(n)$ 。

One thing to note is that when the degree of two polynomials is $n,m$ When , the degree of their product polynomial is $n+m$ , so when calculating using point value representation, you should choose to use the point value representation when calculating the point value representation of two polynomials $n+m+1$ variables, so that the result obtained uniquely represents the product polynomial $C(x)$ 。

However, for a polynomial $A(x)$ For example, substituting any selected variable $(x_0,x_1,x_2,dots,x_{n-1})$ The time complexity of calculating its point value representation is still $O(n^2)$ , and does not have an optimization effect, but the fast Fourier transform solves this problem, reducing the time complexity of converting the coefficient representation into a point value representation as $O(nlogn)$ 。

Fast Fourier Transform FFT

FFT is selected when converting the coefficient representation of the calculation polynomial into the point value representation.complex planeon the unit circleunit complex root $mega_n^k$ Calculate the point value of a polynomial as a variable, here the unit root $mega_n^k$ Satisfy some of the following properties (if you don’t understand anything, you can look up some relevant knowledge about complex numbers by yourself):

$ASIC Flow$

Figure 1 FFT

All of the above properties can be derived fromEuler formulaObtained, the derivation process is not the focus of FFT and is omitted here.

For a polynomial $A(x)=a_0+a_1x+a_2x^2+dots +a_{n-1}x_{n-1}$ , we can divide it to separate even-order terms from odd-order terms, assuming here $n-1$ is an odd number, we get:
$A(x)=(a_0+a_2x^2+dots +a_{n-2}x{n-2})+x(a_1+a_3x2+dots +a_{n-1}x^{n-2})$
We define two polynomials respectively $A_1(x),A_2(x)$ ：
$egin{aligned} &A_1(x)=a_0+a_2x^1+dots +a_{n-2}x^{rac{n}{2}-1} &A_2(x)=a_1+a_3x^2+dots +a_{n-1}x^{rac{n}{2}-1} nd{aligned}$
Then the original polynomial $A(x)$ It can be expressed as:
$A(x)=A_1(x2)+xA_2(x2)$
will $mega_n^k$ Substitute into the above equation:
$egin{aligned} A(mega_nk)&=A_1(\omega_n{2k})+mega_nkA_2(\omega_n{2k}) &=A_1(mega_{rac{n}{2}}k)+\omega_nkA_2(mega_{rac{n}{2}}^k) nd{aligned}$
will $mega_n^{k+rac{n}{2}}$ Substitute into the above equation:
$egin{aligned} A(mega_n{k+\frac{n}{2}})&=A_1(\omega_{n}{2k+n})+mega_n{k+\frac{n}{2}}A_2(\omega_n{2k+n}) &=A_1(mega_n^nimes mega_n{2k})-\omega_nkA_2(mega_n^nimes mega_n^{2k}) &=A_1(mega_n{2k})-\omega_nkA_2(mega_n^{2k}) &=A_1(mega_{rac{n}{2}}k)-\omega_nkA_2(mega_{rac{n}{2}}^k) nd{aligned}$
From this it can be found that $A(mega_n^k)$ and $A(mega_n^{k+rac{n}{2}} )$ Only one symbol is different during the calculation, so when performing enumeration calculations $A(mega_n^k)$ can be obtained directly $A(mega_n^{k+rac{n}{2}} )$ value, using this method to divide and conquer can reduce the complexity to $O(nlogn)$ 。

Inverse Fourier Transform IFFT

Using the above method, the point value representation of the product polynomial is obtained, so the problem that needs to be solved now is how to convert the point value representation back to the coefficient representation.

Assume that the FFT point value of the polynomial is expressed as $(y_0,y_1,y_2,dots ,y_{n-1})$ , its coefficient is expressed as $(a_0,a_1,a_2,dots ,a_{n-1})$ , according to the FFT principle, $y_k$ It can be expressed as follows:
$y_k=um_{i=0}{n-1}a_i(\omega_nk)^i$
take $mega_n^k$ ofconjugate complex numbers $mega_n^{-k}$ , define the vector as follows $(c_0,c_1,c_2,dots ,c_{n-1})$ ：
$c_k=um_{i=0}{n-1}y_i(\omega_n{-k})^i$
Then the following formula can be derived from the definition:

$ASIC Flow$

Figure 2 FFT

For unit roots in the complex plane $mega_n^k$ , has the following properties:
$egin{aligned} um_{i=0}{n-1}(\omega_n{j-k})^i&=0uad(jeq k) um_{i=0}{n-1}(\omega_n{j-k})^i& =mega_n^ 0=1 uad (j= k ) nd{aligned}$
So we can get:
$egin{aligned} &c_k=na_k &a_k=rac{c_k}{n} nd{aligned}$
Therefore, after using FFT to obtain the point value representation of the polynomial, you only need to replace the variable with the conjugate complex number of the originally selected unit root and perform another FFT to obtain the coefficient representation of the polynomial.

There is an issue that needs to be explained here. Our above discussions are all based on** $n$ for $2$ to the power ofunder the conditions, then when $n$ No $2$ Required when the power ofwill $n$ expand to be greater than $n$ the smallest $2$ **, when performing the inverse Fourier transform, we can find through the above derivation that when we calculate $c_k$ in $k$ The value is greater than the original $n$ , does not exist $j=k$ situation, so we obtain $a_k$ for $0$ , representing the polynomial $k$ The coefficient of the secondary term is $0$ 。

But the problem doesn't end here

When the range of some cancer data is very large, when using recursion for calculation, a large amount of recursion will cause stack overflow. Is there any way to avoid recursion?

Iterative implementation of FFT

For such a sequence $(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)$ , we observe the process of dividing it into two parts:

$ASIC Flow$

Figure 3 FFT iteration

We discovered a magical property. After bisecting this sequence, the binary of the sequence can be obtained by flipping the binary of the original sequence. Then we can take advantage of this property and use a $O(n)$ The method can directly obtain the final sequence, thus eliminating the recursive process, and can be implemented by reverse recursion using the final sequence.

Rader algorithm

Rader algorithmThat is an algorithm to implement the above operation, for $N$ number, we take increasing natural numbers $(0,1,2,3,dots)$ It is called a sequential sequence; for each number in the sequential sequence, convert its binary in reverse order and then convert it into decimal, which is called a reverse sequence.
For a sequential sequence, the $i$ The binary number can be regarded as the $i/2$ (Here is integer divisibility) The binary number of the number is shifted one bit to the left, and then according to $i$ The parity of $1$ Or don’t add $1$ 。
Then to get its inverse sequence, you only need to reverse this operation, that is, the $i$ The binary number can be regarded as the $i/2$ The binary number of the number is shifted one bit to the right, and then according to $i$ The parity of $1$ Or don’t add $1$ , where the highest position is the $og_{2}n$ Bit.

Iteratively perform FFT (Butterfly Transform)

The recursive sequence is obtained using Rader algorithm, so how to get the final answer through iteration?
This is actually similar to the idea of iteratively implementing the 01 backpack. $A(x)$ in $n$ When the point value of each power of the unit root is $m=n/2$ The powers of subunit roots are $A_1$ or $A_2$ The point value at has been calculated and stored in $A$ array, then use it directly in the iteration process of the next layer $A$ The answer stored in the array can continue to be iteratively calculated.

Iterative optimization of FFT code implementation

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath> using namespace std; const int N = 300010;
const double PI = acos(-1); int n, m;
struct Complex
{ double x, y; Complex operator+ (const Complex& t) const { return {x + t.x, y + t.y}; } Complex operator- (const Complex& t) const { return {x - t.x, y - t.y}; } Complex operator* (const Complex& t) const { return {x * t.x - y * t.y, x * t.y + y * t.x}; }
}a[N], b[N];
int rev[N], bit, tot; void fft(Complex a[], int inv)
{ for (int i = 0; i < tot; i ++ ) if (i < rev[i]) swap(a[i], a[rev[i]]); for (int mid = 1; mid < tot; mid <<= 1) { auto w1 = Complex({cos(PI / mid), inv * sin(PI / mid)}); for (int i = 0; i < tot; i += mid * 2) { auto wk = Complex({1, 0}); for (int j = 0; j < mid; j ++, wk = wk * w1) { auto x = a[i + j], y = wk * a[i + j + mid]; a[i + j] = x + y, a[i + j + mid] = x - y; } } }
} int main()
{ scanf("%d%d", &n, &m); for (int i = 0; i <= n; i ++ ) scanf("%lf", &a[i].x); for (int i = 0; i <= m; i ++ ) scanf("%lf", &b[i].x); while ((1 << bit) < n + m + 1) bit ++; tot = 1 << bit; for (int i = 0; i < tot; i ++ ) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1)); fft(a, 1), fft(b, 1); for (int i = 0; i < tot; i ++ ) a[i] = a[i] * b[i]; fft(a, -1); for (int i = 0; i <= n + m; i ++ ) printf("%d ", (int)(a[i].x / tot + 0.5)); return 0;
}

FFT (fast Fourier transform)

# Input format

# Output format

# Input and output samples

# Input sample

# Output sample

# Problem analysis

# Introduction to FFT

# Fast Fourier Transform FFT

Figure 1 FFT

# Inverse Fourier Transform IFFT

Figure 2 FFT

# But the problem doesn't end here

# Iterative implementation of FFT

Figure 3 FFT iteration

# Rader algorithm

# Iteratively perform FFT (Butterfly Transform)

# Iterative optimization of FFT code implementation

Input format

Output format

Input and output samples

Input sample

Output sample

Problem analysis

Introduction to FFT

Fast Fourier Transform FFT

Inverse Fourier Transform IFFT

But the problem doesn't end here

Iterative implementation of FFT

Rader algorithm

Iteratively perform FFT (Butterfly Transform)

Iterative optimization of FFT code implementation