做公司网站利润,j建设网站,网站代运营合同,网站小图标素材下载//写在前面 单就FFT算法来说的话#xff0c;下面只给出个人认为比较重要的推导#xff0c;详细的介绍可参考 FFT算法学习笔记 令v[n]是长度为2N的实序列#xff0c;V[k]表示该实序列的2N点DFT。定义两个长度为N的实序列g[n]和h[n]为 g[n]v[2n], h[n]v[2n1], 0n…//写在前面 单就FFT算法来说的话下面只给出个人认为比较重要的推导详细的介绍可参考 FFT算法学习笔记 令v[n]是长度为2N的实序列V[k]表示该实序列的2N点DFT。定义两个长度为N的实序列g[n]和h[n]为 g[n]v[2n], h[n]v[2n1], 0nN 则可进行如下推导 这里所用的FFT算法能够实现O(nlogn)复杂度的离散傅里叶变换和上面最后所得的关系密切相关。 下面进入正题——模意义下的FFT 还是需要先了解一下关于 复序列的DFT的对称性质及一些补充定义 由此可以试想假设说要模的素数p为1e8级别大小那么我们可以把原始的实序列x[n]“拆”一下。 下面假设我们要求的是x[n]卷积y[n]的结果t[n]。 假设q是sqrt(p)级别的一个数我们可以把x[n]/q存到复序列x1[n]的实部x[n]%q存到复序列x1[n]的虚部。这时对x1[n]、y1[n]求DFT再由X1[k]*Y1[k]得到T1[k],整个运算过程中能够产生的最大浮点数为N*q^2级别一般来说还是在可以接受的范围内的。 接下来考虑从卷积结果{T1[k]}中恢复出原始的t[n]的过程。 看一下T1[k]的组成 到这里差不多就可以结束了。发现上面最后一行等号右边有四个“乘积”我们可以把上面四个乘积分别单独拿出来求IDFT就可以恢复出x/y_re/im卷积的结果之后针对不同的乘积考虑需要在模p意义下乘上q^2、q^1或q^0来进行恢复就可以了。 奉上模板 namespace FFT_MO //前面需要有 mod(1e8~1e9级别),upmo(a,b) 的定义
{const int FFT_MAXN118;const db pi3.14159265358979323846264338327950288L;struct cp{db a,b;cp(double a_0,double b_0){aa_,bb_;}cp operator (const cprhs)const{return cp(arhs.a,brhs.b);}cp operator -(const cprhs)const{return cp(a-rhs.a,b-rhs.b);}cp operator *(const cprhs)const{return cp(a*rhs.a-b*rhs.b,a*rhs.bb*rhs.a);}cp operator !()const{return cp(a,-b);}}nw[FFT_MAXN1],f[FFT_MAXN],g[FFT_MAXN],t[FFT_MAXN]; //a-f,b-g,t~c int bitrev[FFT_MAXN]; void fft_init() //初始化 nw[],bitrev[] {int L0;while((1L)!FFT_MAXN) L;for(int i1;iFFT_MAXN;i) bitrev[i]bitrev[i1]1|((i1)(L-1));for(int i0;iFFT_MAXN;i) nw[i]cp((db)cosl(2*pi/FFT_MAXN*i),(db)sinl(2*pi/FFT_MAXN*i));}// n已保证是2的整数次幂 // flag1:DFT | flag-1: IDFTvoid dft(cp *a,int n,int flag1) {int d0;while((1d)*n!FFT_MAXN) d;for(int i0;in;i) if(i(bitrev[i]d))swap(a[i],a[bitrev[i]d]);for(int l2;ln;l1){int delFFT_MAXN/l*flag; // 决定 wn是在复平面是顺时针还是逆时针变化以及变化间距 for(int i0;in;il){cp *leai,*riai(l1);cp *wflag1? nw:nwFFT_MAXN; // 确定wn的起点 for(int k0;k(l1);k){cp ne*ri * *w;*ri*le-ne,*le*lene;le,ri,wdel;}}}if(flag!1) for(int i0;in;i) a[i].a/n,a[i].b/n;}// convo(a,n,b,m,c) a[0..n]*b[0..m] - c[0..nm]void convo(LL *a,int n,LL *b,int m,LL *c){for(int i0;inm;i) c[i]0;int N2;while(Nnm) N1; // N是c扩展后的长度 for(int i0;iN;i) //扩展 a[],b[] 存入f[],g[]注意取模 {LL aain?a[i]:0,bbim? b[i]:0; aa%mod,bb%mod;f[i]cp(db(aa15),db(aa32767));g[i]cp(db(bb15),db(bb32767));}dft(f,N),dft(g,N);for(int i0;iN;i) // 恢复虚部两个“乘积”乘积具体意义见上文{int ji? N-i:0;t[i]((f[i]!f[j])*(!g[j]-g[i])(!f[j]-f[i])*(g[i]!g[j]))*cp(0,0.25);}dft(t,N,-1);for(int i0;inm;i) upmo(c[i],(LL(t[i].a0.5))%mod15);for(int i0;iN;i) // 恢复实部两个“乘积”{int ji? N-i:0;t[i](!f[j]-f[i])*(!g[j]-g[i])*cp(-0.25,0)cp(0,0.25)*(f[i]!f[j])*(g[i]!g[j]);}dft(t,N,-1);for(int i0;inm;i) upmo(c[i],LL(t[i].a0.5)(LL(t[i].b0.5)%mod30));}
} 模板 举个栗子~ hdu 6088 Rikka with Rock-paper-scissors (2017 多校第五场 1004 【组合数学 数论 模意义下的FFT】 //本博客主要参考资料数字信号处理——基于计算机的方法第四版 [美] Sanjit K. Mitra 著 余翔宇 译 转载请注明出处 http://www.cnblogs.com/Just--Do--It/p/7892254.html 谢谢阅读 转载于:https://www.cnblogs.com/Just--Do--It/p/7892254.html
