Fast Fourier transform (FFT)

polynomial

Coefficient representation

If \ (f(x) \) is a \ (n-1 \) degree polynomial, then \ (f(x)=\sum\limits_{i=0}^{n-1}a_i*x^i\)

Where \ (a_i \) is the coefficient of \ (f(x) \), and the complexity of this method is \ (O(n^2) \)

Point value representation

According to the knowledge of elementary school, a polynomial of degree (n-1) can be uniquely determined by points

In other words, if we know the points of a polynomial ((x)_ 1,y_ 1),(x_ 2,y_ 2)…… (x_n,y_n)\)

Then, this polynomial satisfies the unique condition of (y) for any (1) le I / Le n \_ i=\sum\limits_ {j=0}^{n-1}a_ j*x_ i^j\)

So what is the complexity of polynomial multiplication with point values?

First of all, we need to select \ (n \) points, and each point needs to find its value in the polynomial, and the complexity is \ (O(n^2) \)

Then multiply the polynomials represented by the two point values, because \ (c(x_i)=a(x_i)*b(x_i) \), the complexity is \ (O(n) \)

Finally, the interpolation method uses the point value to find the coefficient with the complexity of \ (O(n^2) \) (I can't interpolate yet)

If the point value conversion coefficient and coefficient conversion point value can be realized quickly, wouldn't it be possible to calculate the polynomial multiplication quickly

Pre cheese

complex

If the imaginary number unit \ (I = \ \ sqrt {- 1} \) is defined and \ (a,b \) is a real number, then a number in the form of \ (a+bi \) is called a complex number

Where \ (a \) is the real part of the complex number and \ (b \) is the virtual part of the complex number

In the complex plane, complex number can be expressed as vector, so it has many similar properties with vector. We can also use vector to understand complex number, but complex number has more properties, such as substituting as a number into polynomial

The module length is defined as \ (\ \ sqrt{a^2+b^2} \), and the amplitude is defined as the directed angle from the positive half axis of the axis \ (x \) to the rotation angle of the vector

Complex algorithm:

The addition and subtraction method is the same as vector, with emphasis on multiplication

Geometric definition is: module length multiplication, argument angle addition

Algebra is defined as:

\[(a+bi)*(c+di) \]

\[=ac+adi+bci+bdi^2 \]

\[=(ac-bd)+(ad+bc)i \]

Unit root

We first define the center of the circle as the coordinate origin, and the circle with a radius of \ (1 \) as the unit circle

We divide the circle equally to get the points on the circle. If we take the abscissa of the points on the circle as the real part and the ordinate as the imaginary part, we get the complex numbers

Pictures of online pickpockets (QwQ)

First of all, we don't ask for trouble. Take \ ((1,0) \) as the starting point of \ (n \) points, and record it as \ ({\ Omega_ n} ^ 0 \), counter clockwise \ (K \) is marked as \ ({omega)_ n}^k\)

According to the multiplication of module length and the addition of argument, we can see that_ N ^ k \) is a_ The \ (K \) power of \ ^ 0 \, where \ (\ Omega_ N ^ 1 \) is called the root of the (n \) subunits

According to the argument, we can calculate the value of the_ N ^ k) represents a complex number of \ (cos{\frac{k}{n}2\pi}+i*sin{\frac{k}{n}2\pi} \)

Unit roots have some properties:

\(1.\omega _n^k=cos{\frac{k}{n}2\pi}+i*sin{\frac{k}{n}2\pi}=e^{\frac{2\pi ki}{n}}\)

It is proved that the first step to the second step is obtained by definition, and the second step to the third step is obtained by Euler formula

\(2.\omega _{2n}^{2k}=\omega_{n}^{k}\)

Proof: Omega_ {2n}^{2k}=e^{\frac{2\pi 2ki}{2n}}=e^{\frac{2\pi ki}{n}}=\omega_ {n}^{k}\)

\(3.\omega _{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}\)

Proof: Omega_ {n}^{k+\frac{n}{2}}=\omega _ {2n}^{2k+n}=\omega _ {2n}^{2k}*\omega _ {2n}^{n}=\omega_ {n}^{k}*(cos\pi+i*sin\pi)=-\omega_ {n}^{k}\)

\(4.\omega _{n}^{0}=\omega _{n}^{n}=1\)

Proof: no need to prove it

Before text

This may help you to understand the algorithm:

Fourier, the great immortal, has never seen what a computer looks like. Therefore, the Fourier transform and inverse transformation proposed by him are only methods to transform coefficients into point values and point values into coefficients, and have no effect on reducing the complexity. As for the fast Fourier transform, which is an acceleration method discovered by later generations after studying Fourier transform, it is the optimization of (DFT) and (IDFT)

Discrete Fourier transform (DFT)

Suppose \ (f(x)=\sum\limits_{i=0}^{n-1}a_i*x_i\)

\(DFT(a)=(f(1),f(\omega _{n}),f(\omega _{n}^{2}),......,f(\omega _{n}^{n-1}))\)

To put it more simply, for a polynomial (f(x)), the (1, Omega_ {n},\omega _ {n}^{2},…… ,\omega _ {n} (n-1}) \) to find the point value representation of the polynomial

Inverse discrete Fourier transform (IDFT)

The point value representation of \ (f(x) \) at the root of \ (n \) sub units is transformed into coefficient representation

Here we can answer why we should substitute the root of \ (n \) subunits as \ (x \) polynomials

Suppose (y_0,y_1,y_2,…… y_{n-1}) \) is a polynomial \ (a (x) = \ sum \ limits_ {i=0}^{n-1}a_ i*x_ Discrete Fourier transform of I \)

We have another polynomial \ (b (x) = \ \ sum \ \ limits_ {i=0}^{n-1}=y_ i*x_ I\)

Take the reciprocal of the root of the above-mentioned unit (n) to (1), Omega_ {n}^{-1},\omega _ {n}^{-2},…… ,\omega _ {n} A new discrete Fourier transform (z) is obtained by substituting {- (n-1)}) \) into (B(x) \)_ 0,z_ 1,z_ 2,…… z_{n-1})\)

Then we found that

\[z_k=\sum\limits_{i=0}^{n-1}y_i*(\omega _n^{-k})^i \]

\[=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j*(\omega _n^{i})^j)*(\omega _n^{-k})^i \]

\[=\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(\omega _n^{j-k})^i) \]

For_ {i=0}^{n-1}(\omega _ N ^ {J-K}) ^ I \) let's consider separately:

When \ (j-k=0 \)

The answer is \ (n \)

When \ (j\ne k \)

To sum a series of equal proportions and get \ (\ frac {(\ Omega_ n^{j-k})^n-1}{\omega _ n^{j-k}-1}=\frac{(\omega _ n^n)^{j-k}-1}{\omega _ n^{j-k}-1}=\frac{1^{j-k}-1}{\omega _ n^{j-k}-1}=0\)

therefore

\[\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(\omega _n^{j-k})^i)=n*a_k \]

Namely

\[a_k=\frac{z_k}{n} \]

It is concluded that for the polynomial (B(x)), which takes the discrete Fourier transform of \ (A(x) \) as coefficient, the reciprocal of unit root \ ((1), \ Omega_ {n}^{-1},\omega _ {n}^{-2},…… ,\omega _ {n} Replace ^ {- (n-1)}) \) as \ (x \), and divide the result by \ (n \), which is the coefficient of \ (A(x) \)

This conclusion realizes the transformation of polynomial point values into coefficients

fast fourier transform

We analyzed it and found out the complexity Still \ (O(n^2) \)

The point of learning this crap is to reduce the complexity, so let's talk about how to reduce the complexity

Let's first set \ (a (x) = \ \ sum / limits_ {i=0}^{n-1}a_ i*x_ I\)

We classify the subscripts of (A(x) \) by parity

\[A(x)=(a_0+a_2*x^2+......+a_{n-2}*x^{x-2})+(a_1*x+a_3*x^3+......+a_{n-1}*x^{n-1}) \]

Let's have two polynomials

\[A_1(x)=a_0+a_2*x+......+a_{n-2}*x^{\frac{n}{2}-1} \]

\[A_2(x)=a_1+a_3*x+......+a_{n-1}*x^{\frac{n}{2}-1} \]

Then we can find that \ (A(x)=A_1(x^2)+xA_2(x^2)\)

Set \ (x = \ Omega_ N ^ {K} (k < / frac {n} {2}) \) is replaced by

\[A(\omega _n^k)=A_1(\omega _n^{2k})+\omega _n^kA_2(\omega _n^{2k}) \]

\[=A_1(\omega _{\frac{n}{2}}^{k})+\omega _n^kA_2(\omega _{\frac{n}{2}}^{k}) \]

Set \ (x = \ Omega_ N ^ {K + \ frac {n} {2}} (k < / frac {n} {2}) \) is replaced by

\[A(\omega _n^{k+\frac{n}{2}})=A_1(\omega _n^{2k+n})+\omega _n^{k+\frac{n}{2}}A_2(\omega _n^{2k+n}) \]

\[=A_1(\omega _n^{2k}*\omega _n^n)-\omega _n^kA_2(\omega _n^{2k}*\omega _n^n) \]

\[=A_1(\omega _{\frac{n}{2}}^{k})-\omega _n^kA_2(\omega _{\frac{n}{2}}^{k}) \]

We are not surprised to find that as long as we find \ (A_1(x) \) and (a)_ 2 (x) \) in_ {\frac{n}{2}}^0,\omega _ {\frac{n}{2}}^1,…… ,\omega _ The point value of {frac {n} {2}} ^ {frac {n} {2} - 1}) \) can be used to calculate the value of \ (A(x) \) in ((1, omega)_ {n},\omega _ {n}^{2},…… ,\omega _{n}^{n-1})\)

So we can recursively implement (O(nlogn) \) to solve polynomial multiplication

Note: if we assume \ (f*g=h \), then for both \ (f \) and \ (g \), we need to directly find out the value of \ (2^k \) which is larger than \ (n+m+1 \) (due to the requirement of divide and conquer, the number of points must be the integral power of \ (2 \)

\(code\)

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define mid ((l+r)>>1)
#define eps (1e-8)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=5e6+10;
	const double pi=acos(-1.0);
	int n,m;
	struct complex
	{
		double x,y;
		complex(double tx=0,double ty=0){x=tx,y=ty;}
		inline complex operator + (const complex t) const
		{
			return complex(x+t.x,y+t.y);	
		}
		inline complex operator - (const complex t) const
		{
			return complex(x-t.x,y-t.y);
		}
		inline complex operator * (const complex t) const
		{
			return complex(x*t.x-y*t.y,x*t.y+y*t.x);
		}
	}a[N],b[N];
	inline void fft(int limit,complex *a,int inv)
	{
		if(limit==1) return;
		complex a1[limit>>1],a2[limit>>1];
		for(int i=0;i<limit;i+=2)
		{
			a1[i>>1]=a[i],a2[i>>1]=a[i+1];
		}
		fft(limit>>1,a1,inv);
		fft(limit>>1,a2,inv);
		complex Wn=complex(cos(2.0*pi/limit),inv*sin(2.0*pi/limit)),w=complex(1,0);
		for(int i=0;i<(limit>>1);++i,w=w*Wn)
		{
			a[i]=a1[i]+w*a2[i];
			a[i+(limit>>1)]=a1[i]-w*a2[i];
		}
	}
	inline void main()
	{
		n=read(),m=read();
		for(int i=0;i<=n;++i) a[i].x=read();
		for(int i=0;i<=m;++i) b[i].x=read();
		int limit=1;
		while(limit<=n+m) limit<<=1;
		fft(limit,a,1);
		fft(limit,b,1);
		for(int i=0;i<=limit;++i)
		{
			a[i]=a[i]*b[i];
		}
		fft(limit,a,-1);
		for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].x/limit+0.5));
	}
}
signed main()
{
	red::main();
	return 0;
}

However, we found that it seemed a little slow

Iterative optimization

As we all know, recursion is slow. Is there any way to replace recursion with iteration?

Time for drawing

We can't find that the position of each number after division and division is the binary flip of its position

This law also has a nice name, Butterfly Theorem

So we just need to preprocess the position of each number in the last recursion, and then merge them from bottom to top. Can't we get rid of recursion

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define eps (1e-8)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=5e6+10;
	const double pi=acos(-1.0);
	int n,m;
	int limit=1,len;
	int pos[N];
	struct complex
	{
		double x,y;
		complex(double tx=0,double ty=0){x=tx,y=ty;}
		inline complex operator + (const complex t) const
		{
			return complex(x+t.x,y+t.y);	
		}
		inline complex operator - (const complex t) const
		{
			return complex(x-t.x,y-t.y);
		}
		inline complex operator * (const complex t) const
		{
			return complex(x*t.x-y*t.y,x*t.y+y*t.x);
		}
	}a[N],b[N],buf[N];
	inline void fft(complex *a,int inv)
	{
		for(int i=0;i<limit;++i)
			if(i<pos[i]) swap(a[i],a[pos[i]]);
		for(int mid=1;mid<limit;mid<<=1)
		{
			complex Wn(cos(pi/mid),inv*sin(pi/mid));
			for(int r=mid<<1,j=0;j<limit;j+=r)
			{
				complex w(1,0);
				for(int k=0;k<mid;++k,w=w*Wn)
				{
					buf[j+k]=a[j+k]+w*a[j+k+mid];
					buf[j+k+mid]=a[j+k]-w*a[j+k+mid];
				}
			}
			for(int i=0;i<limit;++i) a[i]=buf[i];
		}
	}
	inline void main()
	{
		n=read(),m=read();
		for(int i=0;i<=n;++i) a[i].x=read();
		for(int i=0;i<=m;++i) b[i].x=read();
		while(limit<=n+m) limit<<=1,++len;
		for(int i=0;i<limit;++i)
			pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
		fft(a,1);
		fft(b,1);
		for(int i=0;i<=limit;++i) a[i]=a[i]*b[i];
		fft(a,-1);
		for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].x/limit+0.5));
	}
}
signed main()
{
	red::main();
	return 0;
}

Butterfly operation

Think about it

for(int r=mid<<1,j=0;j<limit;j+=r)
{
	complex w(1,0);
	for(int k=0;k<mid;++k,w=w*Wn)
	{
		buf[j+k]=a[j+k]+w*a[j+k+mid];
		buf[j+k+mid]=a[j+k]-w*a[j+k+mid];
	}
}
for(int i=0;i<limit;++i) a[i]=buf[i];

The reason for adding the array of \ (buf \) is that the value of \ (a \) assigned twice will change. We can store the value of array \ (a \) in advance, and then optimize the array of \ (buf \)

for(int k=0;k<mid;++k,w=w*Wn)
{
	complex x=a[j+k],y=w*a[j+k+mid];
	a[j+k]=x+y;
	a[j+k+mid]=x-y;
}

Three times to two times optimization

Looking at the code above, we ran the fat rabbit three times. Now we have a way to run one less time

Suppose we find \ (f(x)*g(x) \)

Let complex polynomial \ (h(x)=f(x)+i*g(x) \, real part \ (f(x) \, virtual part \ (g(x) \)

Then \ (h(x)^2=(f(x)+i*g(x))^2=f(x)^2-g(x)^2+i*2*f(x)*g(x) \)

We just divide the imaginary part of \ (h(x)^2 \) by \ (2 \)

Full version:

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define eps (1e-8)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=5e6+10;
	const double pi=acos(-1.0);
	int n,m;
	int limit=1,len;
	int pos[N];
	struct complex
	{
		double x,y;
		complex(double tx=0,double ty=0){x=tx,y=ty;}
		inline complex operator + (const complex t) const
		{
			return complex(x+t.x,y+t.y);	
		}
		inline complex operator - (const complex t) const
		{
			return complex(x-t.x,y-t.y);
		}
		inline complex operator * (const complex t) const
		{
			return complex(x*t.x-y*t.y,x*t.y+y*t.x);
		}
	}a[N];
	inline void fft(complex *a,int inv)
	{
		for(int i=0;i<limit;++i)
			if(i<pos[i]) swap(a[i],a[pos[i]]);
		for(int mid=1;mid<limit;mid<<=1)
		{
			complex Wn(cos(pi/mid),inv*sin(pi/mid));
			for(int r=mid<<1,j=0;j<limit;j+=r)
			{
				complex w(1,0);
				for(int k=0;k<mid;++k,w=w*Wn)
				{
					complex x=a[j+k],y=w*a[j+k+mid];
					a[j+k]=x+y;
					a[j+k+mid]=x-y;
				}
			}
		}
	}
	inline void main()
	{
		n=read(),m=read();
		for(int i=0;i<=n;++i) a[i].x=read();
		for(int i=0;i<=m;++i) a[i].y=read();
		while(limit<=n+m) limit<<=1,++len;
		for(int i=0;i<limit;++i)
			pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
		fft(a,1);
		for(int i=0;i<=limit;++i) a[i]=a[i]*a[i];
		fft(a,-1);
		for(int i=0;i<=n+m;++i) printf("%d ",(int)(a[i].y/limit/2+0.5));
	}
}
signed main()
{
	red::main();
	return 0;
}

Please consider whether to use it according to the range of the subject value

Reference blog
%%%attack
%%%rabbithu

Some examples

Irregular update
Logue P1919 A*B Problem upgrade
Luogu P3338 [ZJOI2014] force

Tags: less

Posted on Mon, 29 Jun 2020 04:04:11 -0400 by Wade