【BZOJ2194】快速傅里叶之二，FFT和一点奇怪的想法

传送门(权限题)
题意：给定数列a,b,求ck=∑i=kn−1aibi−k\displaystyle c_k=\sum_{i=k}^{n-1}a_{i}b_{i-k},n≤105n\leq10^5
思路：
（偶然发现一点新东西的奇妙感，如果神牛们早就知道了就无视掉我这篇胡扯的博文吧）
下面的n默认为是2的正整数次幂
传统FFT模板题，网上的传统做法是把a或b翻转过来，化成卷积形式然后直接做就可以了
但是我在做这道题的时候并没有想到这种方法，而是一直在考虑在多项式方面的问题，从一般感觉上来说，−(i−k)+i=k-(i-k)+i=k，那怎么化成负的呢？
我们知道一般形式为A(x)=a0+a1x+a2x2+...+an−1xn−1A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}，称其为x的一个次数界为n的多项式；对于多项式乘法，设结果为c，相乘的两个多项式为a,b;ai,bi,cia_i,b_i,c_i分别表示a,b,c次数为i的项，那么ci=∑j=0iajbi−jc_i=\displaystyle \sum_{j=0}^i a_jb_{i-j}，这是一个卷积形式，这也是FFT是处理卷积问题的利器的原因。回到原来的问题上，我们从多项式的项次数上考虑，那么我们要求的ckc_k不就是”a的次数为i的项乘b的次数为-(i-k)的项”的累和？
抱着试一试与乱搞的心态，我写了一个这样的东西
A(x)=a0+a1x+a2x2+...+an−1xn−1A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}
B(x)=b0+b1x1+b2x2+...+bn−1xn−1B(x)=b_0+b_1x^1+b_2x^2+...+b_{n-1}x^{n-1}
C(x)=A(x)B(x−1)C(x)=A(x)B(x^{-1})
那么多项式C(x)C(x)的0~n-1次项的系数就是答案
处理B(x)B(x)的FFT的时候，像逆DFT一样，把单位复数根ωn\omega_n变成ω−1n\omega_n^{-1}就好了
但是C(x)C(x)里有负指数项，正确答案在哪里呢？
抱着弃疗的心态输出了c0,c1..cn−1c_0,c_1..c_{n-1}，发现样例过了，对拍过了，然后就A了
静下心来仔细思考，我们在算出A,B的DFT，再做完向量乘法时，实际上的C是这样的
ci=∑0<j<n,0<j+i<najbj+i\displaystyle c_i=\sum_{0
C(x)=c−(n−1)x−(n−1)+c−(n−2)x−(n−2)+..+c0+c1x+c2x2+...cn−1xn−1C(x)=c_{-(n-1)}x^{-(n-1)}+c_{-(n-2)}x^{-(n-2)}+..+c_0+c_1x+c_2x^2+...c_{n-1}x^{n-1}
但在做逆DFT时，因为逆DFT是插值法，求出来的形式一定与多项式的一般形式相符合而不会出现负指数，而插值正好是2n次单位复数根的k(k=0~2n-1)次方，把负号移到它们头上，所求出来的项的次数是在mod 2n\mod\ 2n意义下的，也就是这些负指数项正好加上2n(我们做FFT时默认A,B,C次数为2n，不过它们实际最高次数并不一定是2n)
上面可能说的很模糊，举例来说，我们在C(x)C(x)的原式中所谓的第−j-j项c−jx−jc_{-j}x^{-j}，因为x=ω02n,ω12n..ω2n−12nx=\omega_{2n}^0,\omega_{2n}^1..\omega_{2n}^{2n-1}，所以这里的c−jx−j=c−jx2n−jc_{-j}x^{-j}=c_{-j}x^{2n-j},显然2n−j>02n-j>0
说白了，我们对A(x)B(x−1)A(x)B(x^{-1})的2n次单位复数根处插值时，得到的多项式实际是
C′(x)=c0+c1x+c2x2+...+cn−1xn−1+c−(n−1)xn+1+...+c−2x2n−2+c−1x2n−1C'(x)=c_0+c_1x+c_2x^2+...+c_{n-1}x^{n-1}+c_{-(n-1)}x^{n+1}+...+c_{-2}x^{2n-2}+c_{-1}x^{2n-1}
不得不说，这样做的基础就是n次单位复数根在乘法下构成n阶群以及插值法确定多项式的唯一性
（上面这句话是我瞎扯淡的，请无视）
代码：

#include<cstdio>
#include<complex>
using namespace std;
int n;
int pos[1<<18],bit[1<<18];
const double pi=acos(-1.0);
typedef complex<double> Co;
Co a[1<<18],b[1<<18];
void FFT(Co a[],bool tp,int len)
{for (int i=0;i<len;++i)if (pos[i]>i) swap(a[pos[i]],a[i]);for (int i=2;i<=len;i=i<<1){Co w(cos(2*pi/i),(tp?-1:1)*sin(2*pi/i)),wn,p,q;for (int j=0;j<len;j+=i){wn=Co(1,0);for (int k=0;k<i/2;++k)p=a[k+j],q=wn*a[k+j+i/2],a[k+j]=p+q,a[k+j+i/2]=p-q,wn=w*wn;}}
}
main()
{scanf("%d",&n);for (int x,y,i=0;i<n;++i)scanf("%d%d",&x,&y),a[i]=Co(x,0),b[i]=Co(y,0); int wei=log2(n*2)+1,len=1<<wei;for (int t=1,i=0;i<wei;++i,t+=t) bit[t]=i;for (int i=1;i<len;++i) pos[i]=pos[i-(i&-i)]|(1<<wei-bit[i&-i]-1);FFT(a,0,len);FFT(b,1,len);for (int i=0;i<len;++i) a[i]=a[i]*b[i];FFT(a,1,len);for (int i=0;i<n;++i) printf("%d\n",(int)(a[i].real()/len+0.5));
}

【BZOJ2194】快速傅里叶之二，FFT和一点奇怪的想法相关推荐

bzoj2194 快速傅里叶之二
题意:对于k = 0 ... n求解: 首先把i变成从0开始我们发现a和b的次数(下标)是成正比例的,这不可,于是反转就行了. 反转b的话,会发现次数和是n + k,这不可. 反转a就很吼了. 这 ...
前端必知必会HTTP请求系列（二）简单一点的HTTP协议
http协议用户客户端和服务器之间的通信 http协议和TCP/IP协议族内的其他众多协议相同,用于客户端和服务器之间的通信. 那么问题来个如果两台服务器之间一台服务器向另一台服务器进行接口请求那谁是 ...
基二FFT时间抽取和频域抽取算法
目录: 一.理解离散傅立叶变换(FFT结果的物理意义(编程参考用)) 二.基二频域算法原理三.基二FFT的C语言实现(时域) 四.基二FFT的C语言实现(频域) ------------------ ...
基于STM32F407实现快速傅里叶变化（FFT），计算指定频率的幅值
本人的课题是关于EIT采集系统,简单的说就是往人体注入特定频率的恒流源,再采集电压信号,通过分析电阻抗分布进行成像.采集的电压信号是需要进行FFT处理,只保留注入频率的信号成分.本文主要介绍如何在ST ...
快速傅里叶算法（FFT）快在哪里？
目录前言 1.DFT算法 2.FFT算法 2.1 分类 2.2 以基2 DIT(时间抽取) FFT 算法为例 2.2.1 一次分解 2.2.2 多次分解参考前言对信号分析的过程中,为了能换 ...
正则表达式的一点奇怪
(\d?) 11a111d1d111d1 -- 结果是 15处匹配, -- 可以推断是,用正则原子一个一个的去匹配字符串,然后得出一个结论. 如果正则是(\d) -- 结果是 10处匹配 \s匹配空格 ...
[转载]关于申请国外博后的一点经验和想法
很多人问,申请博后难不难啊.这其实是个小马过河的问题,这取决于你自己到底是大牛还是小松鼠了. 如果要申请博后,我觉得首先你需要给你自己一个定位,定位标准是这样的,分3类:1.你的博士老板是牛老板:2. ...
带三维团队半年的一点总结和想法
文章版权由作者李晓晖和博客园共有,若转载请于明显处标明出处:http://www.cnblogs.com/naaoveGIS/ 1. 一个有些突然的人事变动公司二维团队和三维团队由于历史原因, ...
关于申请国外博后的一点经验和想法
很多人问,申请博后难不难啊.这其实是个小马过河的问题,这取决于你自己到底是大牛还是小松鼠了. 如果要申请博后,我觉得首先你需要给你自己一个定位,定位标准是这样的,分3类:1.你的博士老板是牛老板:2. ...

【BZOJ2194】快速傅里叶之二，FFT和一点奇怪的想法

【BZOJ2194】快速傅里叶之二，FFT和一点奇怪的想法相关推荐

最新文章

热门文章