洛谷p3803 FFT入门

洛谷p3803 FFT入门

浅谈FFT的基本原理：
对于两个多项式的卷积，可以O(n^2)求出来（妥妥的暴力）
显然一个多项式可以用$a_0+a_{1}X+a_2{X}^2+a_3{X}^3+a_4{X}^4$……表示。
也可以用(x1,y1),(x2,y2),(x3,y3),(x4,y4)的点集来表示。
用点值表示有一个好处：两个多项式的卷积可以直接取相同的x值，y值相乘得到。
那么，怎么转化为点值表示呢？
直接代进去？显然也是O(n^2)，没用……
$$\begin{array}{rl}设A(X)& =a_0+a_{1}X+a_2{X}^2+a_3{X}^3+a_4{X}^4+a_5{X}^5..\\ A_0(X)& =a_0+a_{2}X+a_4{X}^2..\\ A_1(X)& =a_1+a_{3}X+a_5{X}^2..\\ & \\ 那么A(X)& =A0(X^2)+X\ast A1(X^2)\end{array}$$

然而，$X^2$仍会有N种不同的取值。
引入一个复数根可以很好的解决这个问题：设$i=\sqrt{-1}$
$$\begin{gathered} \boxed{e^{i\theta }=cos\theta +i\ast sin\theta \\ (cos\theta +i\ast sin\theta {)}^k=cosk\theta +i\ast sink\theta } \end{gathered}$$

$$\begin{array}{rl}记W_n^1& =cos\frac{2\pi }{n}+i\ast sin\frac{2\pi }{n}\\ W_n^k& =(W_n^1{)}^k=cos\frac{2\pi k}{n}+i\ast sin\frac{2\pi k}{n}表示一个单位圆分成n份中的第k份(k次方)\\ (W_n^k{)}^2& =(W_n^1{)}^{2k}=cos\frac{2\pi \ast 2k}{n}+i\ast sin\frac{2\pi \ast 2k}{n}=cos\frac{2\pi k}{\frac{n}{2}}+i\ast sin\frac{2\pi k}{\frac{n}{2}}=W_{\frac{n}{2}}^k\\ & \\ 那么A(W_n^k)& =A_0((W_n^k{)}^2)+W_n^k\ast A_1((W_n^k{)}^2)\\ & =A_0(W_{\frac{n}{2}}^k)+W_n^k\ast A_1(W_{\frac{n}{2}}^k)\end{array}$$

根据三角函数的周期性：2π为一个周期。
$$\begin{array}{rl}{W}_n^k& =W_n^{kmodn}\\ W_n^{\frac{n}{2}}& =W_2=cos\frac{2\pi }{2}+i\ast sin\frac{2\pi }{2}=cos\pi +i\ast sin\pi =-1\\ W_n^{t+\frac{n}{2}}& =W_n^t{W}_n^{\frac{n}{2}}=-W_n^t\\ & \\ A(W_n^k)& =A_0(W_{\frac{n}{2}}^k)+W_n^k\ast A_1(W_{\frac{n}{2}}^k)\\ 可以拆成两部分,设t& <\frac{n}{2}\\ A(W_n^t)& =A_0(W_{\frac{n}{2}}^t)+W_n^t\ast A_1(W_{\frac{n}{2}}^t)\\ A(W_n^{t+\frac{n}{2}})& =A_0(W_{\frac{n}{2}}^{t+\frac{n}{2}})+W_n^{t+\frac{n}{2}}\ast A_1(W_{\frac{n}{2}}^{t+\frac{n}{2}})\\ & =A_0(W_{\frac{n}{2}}^t)+W_n^t{W}_n^{\frac{n}{2}}\ast A_1(W_{\frac{n}{2}}^t)\\ & =A_0(W_{\frac{n}{2}}^t)-W_n^t\ast A_1(W_{\frac{n}{2}}^t)\end{array}$$

$$A(W_n^t)和A(W_n^{t+\frac{n}{2}})可通过A_0(W_{\frac{n}{2}}^t)和W_n^t\ast A_1(W_{\frac{n}{2}}^t)推导得到$$

证毕，啦啦啦！

往下递归 $W_{\frac{n}{2}}$的答案。只需要 $lo{g}_{2}n$层，一直递归下去算就可以将时间优化到$O(nlo{g}_{2}n)$

分治过程图解

非递归的蝴蝶变换

　上图摘自算法导论。算法导论下载，提取码qcfs

void FFT(Complex *A,int n,int DFT)
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(A[i],A[rev[i]]);
    for(int s=1;(1<<s)<=n;s++)
    {
        int mi=(1<<s);
        Complex wn=Complex(cos(2*pi/mi),sin(2*pi/mi)*DFT);
        for(int t=0;t<n;t+=mi)
        {
            Complex w=Complex(1,0);
            for(int j=0;j<(mi>>1);j++)
            {
                Complex u=A[t+j],v=w*A[t+j+(mi>>1)];
                A[t+j]=u+v; A[t+j+(mi>>1)]=u-v;
                w=w*wn;
            }
        }
    }
    if(DFT==-1) for(int i=0;i<n;i++) A[i].real/=n,A[i].image/=n;
}
for(int i=0;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));

一题洛谷的裸题p3803
【代码】

#include<cstdio>
#include<algorithm>
#include<cmath>

using namespace std;

typedef long long ll;

const int N=1e6+10;
const double pi=acos(-1.0);
int n,m,nn,k;
struct Complex {
	double real,image;   //real+i*image
	Complex(){}
	Complex(double _real,double _image)
	{
		real=_real; image=_image;
	}
	friend Complex operator + (Complex A,Complex B) { return (Complex(A.real+B.real,A.image+B.image)); }
	friend Complex operator - (Complex A,Complex B) { return (Complex(A.real-B.real,A.image-B.image)); }
	friend Complex operator * (Complex A,Complex B) { return (Complex(A.real*B.real-A.image*B.image,A.real*B.image+A.image*B.real)); }
} A[N<<2],B[N<<2],C[N<<2];
int rev[N<<2];

void FFT(Complex *A,int n,int DFT)
{
	for(int i=0;i<n;i++) if(i<rev[i]) swap(A[i],A[rev[i]]);
	for(int s=1;(1<<s)<=n;s++)
	{
		int mi=(1<<s);
		Complex wn=Complex(cos(2*pi/mi),sin(2*pi/mi)*DFT);
		for(int t=0;t<n;t+=mi)
		{
			Complex w=Complex(1,0);
			for(int j=0;j<(mi>>1);j++)
			{
				Complex u=A[t+j],v=w*A[t+j+(mi>>1)];
				A[t+j]=u+v; A[t+j+(mi>>1)]=u-v;
				w=w*wn;
			}
		}
	}
	if(DFT==-1) for(int i=0;i<n;i++) A[i].real/=n,A[i].image/=n;
}

int main()
{
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) { int xx; scanf("%d",&xx); A[i]=Complex(xx,0); }
	for(int i=0;i<=m;i++) { int xx; scanf("%d",&xx); B[i]=Complex(xx,0); }
	m+=n; k=0;
	for(nn=1;nn<=m;nn<<=1) k++;
	for(int i=0;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
	FFT(A,nn,1); FFT(B,nn,1);
	for(int i=0;i<=nn;i++) C[i]=A[i]*B[i];
	FFT(C,nn,-1);
	for(int i=0;i<=m;i++) printf("%d ",(int)(C[i].real+0.5));
	printf("\n");
	return 0;
}

ps：这篇总结写得真心累……