洛谷P3803 fft模板

概念：

先引入多项式的两个表示方法：

$1.$系数表示法，唯一确定一个$n$次的多项式需要每一项的系数，从低次项到高次项依次写成向量的形式，这个向量就能唯一确定一个多项式。即向量$r=[a_0,a_1,.....,a_n]$，就能唯一确定多项式$f(x)=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^n$。

此时计算一次$f(x)$的时间复杂度是$O(n)$。多项式相加结果即向量相加，复杂度为$O(n)$，多项式相乘即向量相乘，时间复杂度是$O(n^2)$。

等长度向量相乘的方法：

设向量$x=[a_0,a_1,...,a_n]$与向量$y=[b_0,b_1,...,b_n]$相乘，得到的还是长度为$n$的向量$z=[c_0,c_1,...,c_n]$，$c_k$为$x,y$中所有下标之和为$k$的对的积的和，即$c_k={\sum }_{i=0}^k{a}_i{b}_{k-i}$。

$2.$点值表示法，唯一确定一个$n$次的多项式需要$n+1$个点（可以理解为解$n+1$元方程需要$n+1$个方程，每一个点可以确定所有系数的一组关系，有$n+1$个系数，就是$n+1$元方程）。点是二维的，所以需要一对列表（此处非向量）来表示这个多项式，即$[x_1,x_2,...,x_n]$，与$[y_1,y_2,...,y_n]$就可以唯一确定。此时多项式相加依然满足$y$直接相加，乘法直接相乘即可。此时多项式加减法的时间复杂度都是$O(n)$。

单位根的概念：

满足方程$z^n=1$的复数解$z$即$n$次单位根。

性质$1$：单位根的模一定是$1$，并且$n\theta =2k\pi$

由欧拉定理$e^{i\theta }=cos\theta +isin\theta$可知，任意一个右式形式的复数都可以被写成左式的形式，同时注意到，右式的模$\sqrt{co{s}^2\theta +si{n}^2\theta =1}$，这是一个长度为$1$的向量的形式，更一般的，假设长度为$r$，那么形式应是$r{e}^{i\theta }=rcos\theta +risin\theta$的形式（在复平面，$r$即是向量的模，$\theta$是向量的辐角）。因此，设$z=r{e}^{i\theta }=rcos\theta +risin\theta \Rightarrow z^n=r^n{e}^{in\theta }=r^{n}cosn\theta +r^{n}isinn\theta =1\Rightarrow r^n(cosn\theta +isinn\theta )=1$，右式不存在有虚数，因此$si{n}^{n}n\theta =0\Rightarrow n\theta =2k\pi$，代入有$r^n=1\Rightarrow r=1$。

性质$2：n$次单位根有$n$个，并且第$k$个单位根写作$w_n^k$，有$w_n^k=r{e}^{2\pi \frac{k}{n}}$，单位根是从$0$标号到$n-1$的。

注意到每个单位根与自己相乘就是让自己的辐角翻倍，因此，单位根是逆时针顺序在单位圆上排序的，间隔角度是$2\pi /n$，那么第$k$个就是原命题的形式了。

性质$3$：$w_{dn}^{dk}=w_n^k$，$w_n^{k+\frac{n}{2}}=w_{\frac{n}{2}}^k$

展开计算，上面的性质是容易得到的。

算法：已知两多项式系数，求两多项式的系数向量乘积

已知多项式的系数向量分别为：$x=[a_0,a_1,...,a_n]$与$y=[b_0,b_1,...,b_n]$，求乘积向量，其中$c_i={\sum }_{j=0}^i{a}_j{b}_{i-j}$，时间复杂度是$O(n^2)$，难以接受。注意到点值式的乘积运算时$O(n)$的，那么可以利用点值式来计算，现在的问题是，系数如何转换为点值，也就是该用拿几个点？求出点值合并完成后如何转换为系数？并且计算一次值的时间复杂度是$O(n)$，计算$n$项的时间复杂度同样是$O(n^2)$，该怎么加快计算点值？

上述几个问题能被单位根完美解决，我们选用的点值就是单位根上的点值。

设原来的向量为$A=[a_0,a_1,...,a_n]$

我们把系数按奇偶次序分开两个向量，得到：

$A^{[0]}=[a_0,a_2,a_4....]$

$A^{[1]}=[a_1,a_3,a_5....]$

对应的多项式是：

$f(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+....$

$f^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+....$

$f^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+....$

把两个形式变换得到：

$f(x)=f^{[0]}(x^2)+x{f}^{[1]}(x^2)$

假设我们要计算单位根$w_n^k$

带入得到：

$f(w_n^k)=f^{[0]}((w_n^k{)}^2)+w_n^k{f}^{[1]}((w_n^k{)}^2)$

化简得到：

$f(w_n^k)=f^{[0]}(w_{\frac{n}{2}}^k)+w_n^k{f}^{[1]}(w_{\frac{n}{2}}^k)$

注意到单位根的性质$2$，我们求的如果是$w_n^{k+\frac{n}{2}}$，最后会是

$f(w_n^{k+\frac{n}{2}})=f^{[0]}(w_{\frac{n}{2}}^k)-w_n^k{f}^{[1]}(w_{\frac{n}{2}}^k)$

因此，我们只需要计算出前一半的单位根对应的$f^{[0]}$和$f^{[1]}$，后一半是可以$O(1)$得到的。

同时前式和后式是一样的问题，于是可以分治解决，最后合并。

最后给出点值转换回系数的方法：

$a_i=\frac{c_i}{n}$

由于$fft$使用$double$保存复数向量，所以最后要四舍五入，这里也给出一种四舍五入的方式：$(int)(a+0.5)$

由于算法基于对称的方法，所以只适合规模是$2$的幂次的情形，因此，如果规模不为$2$的幂次，需要添加额外项使规模是$2$的幂，添加的额外项系数是0。

#include<bits/stdc++.h>
#define ll long long
#define complex Complex
using namespace std;

int read()
{
    int ret=0,base=1;
    char ch=getchar();
    while(!isdigit(ch))
    {
        if(ch=='-') base=-1;
        ch=getchar();
    }
    while(isdigit(ch))
    {
        ret=(ret<<3)+(ret<<1)+ch-48;
        ch=getchar();
    }
    return ret*base;
}

const double pi=acos(-1);

struct complex
{
    double a,b;
    complex(){a=0;b=0;}
    complex(double x,double y){a=x;b=y;}
    complex operator+(complex x){return complex(a+x.a,b+x.b);}
    complex operator-(complex x){return complex(a-x.a,b-x.b);}
    complex operator*(complex x){return complex(a*x.a-b*x.b,a*x.b+b*x.a);}
};

void fft(vector<complex>&a,int limit,int op)
{
    if(limit==1) return;
    vector<complex>ret1(limit>>1),ret2(limit>>1);
    for(int i=0;i<limit;i+=2) ret1[i>>1]=a[i],ret2[i>>1]=a[i+1];
    fft(ret1,limit>>1,op);fft(ret2,limit>>1,op);
    complex base=complex(cos(2.0*pi/limit),op*sin(2.0*pi/limit)),now=complex(1,0);
    for(int i=0;i<(limit>>1);i++)
    {
        a[i]=ret1[i]+now*ret2[i];
        a[i+(limit>>1)]=ret1[i]-now*ret2[i];
        now=now*base;
    }
}

int n,m;
vector<complex>a(2100005),b(2100005);

int main()
{
    n=read(),m=read();
    for(int i=0;i<=n;i++) a[i].a=read();
    for(int i=0;i<=m;i++) b[i].a=read();
    int limit=1;
    while(limit<=n+m) limit<<=1;
    fft(a,limit,1);fft(b,limit,1);
    for(int i=0;i<=limit;i++) a[i]=a[i]*b[i];
    fft(a,limit,-1);
    for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].a/limit+0.5));
    return 0;
}