我学科技?真的假的?
起因是在做斯特林数,发现没事就要佛佛塔或者呢塔塔,遂学之。
FFT
现在要求两个多项式的卷积。
卷积
现在给你两个序列 \(a_i,b_i\) 并定义 \(c_i\) 为
\[c_i=\sum_{p+q=i}a_ib_i \]可以认为 \(c_i\) 就是 \(a_i,b_i\) 的卷积。换句话说就是给你两个多项式 \(F(x)=\sum_{i=0}^n a_ix^i,G(x)=\sum_{i=0}^nb_ix^i\),然后 \(H(x)=F(x)*G(x)=\sum_{i=0}^{2*n}c_ix^i\),手玩一下二次式乘二次式,算出来的多项式就是卷积结果,多项式的 \(i\) 次式就是 \(c_i\)。
点值表示法
现在让你算两个多项式的乘法就可以 \(O(n^2)\) 算了,但是不够快,考虑优化这个过程,但是发现难以优化。
前人的智慧指出:对任意 \(n\) 个数点(x不同)都可以唯一地插值出来一个对应的 \(n\) 次多项式,既然 \(H=(F*G)\) 则对于原先的数点必定有 \(H(x_0)=F(x_0)*G(x_0)\)。确立这些 \(H(x)\) 的数点的过程应当是 \(O(n)\) 的,如果现在可以用一个比较高速的方法插值出 \(H(x)\) 就可以实现复杂度优化了。
单位根
如果你不是我且在看这段建议先学复数及其复平面上的三角表示和加减乘除于复平面上的表示。
考虑在复平面内给单位圆劈成 \(n\) 份,\(x\) 正方向向上的第一刀对应的向量 \(\omega_n\) 肯定满足 \(\omega_n^n=1\),因为转了一圈,就把这个东西叫单位根。
这个东西肯定就有一些性质比如
\[\omega_{2n}^{2k}=\omega_n^k \]砍两倍次数每次拼两个上去和砍一倍每次拼一个一样...
\[\omega_{2n}^k=-\omega_{2n}^{k+n} \]因为转了一半过去。
分治法法塔(FFT)
现在有这样一个多项式
\[F(x)=a_0+a_1*x+a_2*x^2+...a_{n-1}*x^{n-1} \]我直接给他奇偶劈两半,不妨认为 \(n\) 是偶数。
\[=(a_0+a_2*x^2+...a_{n-2}*x^{n-2})+(a_1*x+a_3*x^3+...a_{n-1}*x^{n-1}) \]再设
\[F_1(x)=a_0+a_2*x+a_4*x^2+...a_{n-2}*x^{\frac{n}{2}-1} \]\[F_2(x)=a_1+a_3*x+a_5*x^2+...+a_{n-1}*x^{\frac{n}{2}-1} \]然后就有
\[F(x)=F_1(x^2)+xF_2(x^2) \]那很明显这个 \(x\) 太丑了直接给他换成幂次等效于旋转的单位根。
\[F(\omega_n^k)=F_1(\omega_n^{2k})+\omega_n^k F_2(\omega_n^{2k}) \]\[=F_1(\omega_{n/2}^{k})+\omega_n^k F_2(\omega_{n/2}^{k}) \]同理
\[F(\omega_n^{k+n/2})=F_1(\omega_{n/2}^{k})-\omega_n^k F_2(\omega_{n/2}^{k}) \]你知道这个玩意有多牛逼吗就是说 \(F(\omega_n^k)\) 和 \(F(\omega_n^{k+n/2})\) 是可以同时算出来的,加号改减号即可而且这个玩意跨度正好是 \(n/2\)。
直接考虑一个递归分治状物然后根据6定理得到时间复杂度为 \(O(nlogn)\)。
蛋是这个玩意常数大到飞起是无法通过板题的,至少复杂度对了常数优化都是后事。
先考虑这样一个问题:这样kuku算完算出来的是一堆 \(F(x)\) 啊插值的事情还没解决。
挨法法塔(IFFT)
现在要解决如何把点值表示法转换为一般的系数表示法,而且要充分利用单位根优势,要不然不如去考虑拉插求系数。
然后前人的智慧再次发力提出一个牛逼的构造,定义
\[c_k=\sum_{i=0}^{n-1}y_i(\omega_{n}^{-k})^i \]\(y_i\) 就是刚才算出来的那坨 \(F()\)。
然后这个东西可以化。
\[c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i \]\[=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i \]\[=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i \]然后先放在这,设 \(S(x)=\sum_{i=0}^{n-1}x^i\) 并带入 \(\omega_n^k\) 然后构造等比数列求和式
\[\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-1 \]\[S(\omega_n^k)=0 \]但是这是建立在 \(\omega_n^{kn}=1\) 的基础上的,\(k=0\) 时原式为 \(n\)
所以这个
\[\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i \]当且仅当 \(j=k\) 时为 \(n\),否则为 \(0\)。
进而 \(a_k=\frac{c_k}{n}\)。
所以流程就是:FFT,点乘计算,反着FFT(带个-1),然后除以 \(n\),即次数。
放一下板子代码方便以后照搬。
#include<bits/stdc++.h>
#define db double
#define MAXN 10000005
using namespace std;
int n,m,pw,len=1,loc[MAXN];
const db pi=acos(-1.0);
struct comp{
db x,y;
comp(db a=0,db b=0){x=a,y=b;}
}a[MAXN],b[MAXN];
comp operator+(comp a,comp b){return comp(a.x+b.x,a.y+b.y);}
comp operator-(comp a,comp b){return comp(a.x-b.x,a.y-b.y);}
comp operator*(comp a,comp b){return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline void FFT(comp *A,int typ){
for(int i=0;i<len;i++)if(i<loc[i])swap(A[i],A[loc[i]]);
for(int mid=1;mid<len;mid<<=1){
comp wn(cos(pi/mid),typ*sin(pi/mid));
for(int r=mid<<1,j=0;j<len;j+=r){
comp w(1,0);
for(int k=0;k<mid;k++,w=w*wn){
comp x=A[j+k],y=w*A[j+k+mid];
A[j+k]=x+y;
A[j+k+mid]=x-y;
}
}
}
}
signed 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(len<=n+m)len<<=1,++pw;
for(int i=0;i<len;i++)loc[i]=(loc[i>>1]>>1)|((i&1)<<(pw-1));
FFT(a,1);
FFT(b,1);
for(int i=0;i<=len;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/len+0.5));
return 0;
}
正在施工。
标签:多项式,sum,FFT,db,comp,相关,omega From: https://www.cnblogs.com/Claimo0611/p/18677569