FFT
2016-06-01 08:36
225 查看
FFT新手请参考下面链接
http://www.gatevin.moe/acm/fft算法学习笔记/
http://www.gatevin.moe/acm/fft算法学习笔记/
#include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #include<cmath> using namespace std; const int maxn=220005; const double eps(1e-8); typedef long long LL; const double PI = acos(-1.0); struct Complex { double real, image; Complex(double _real, double _image) { real = _real; image = _image; } Complex(){} }; Complex operator + (const Complex &c1, const Complex &c2) { return Complex(c1.real + c2.real, c1.image + c2.image); } Complex operator - (const Complex &c1, const Complex &c2) { return Complex(c1.real - c2.real, c1.image - c2.image); } Complex operator * (const Complex &c1, const Complex &c2) { return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real); } int rev(int id, int len) { int ret = 0; for(int i = 0; (1 << i) < len; i++) { ret <<= 1; if(id & (1 << i)) ret |= 1; } return ret; } Complex A[maxn<<1]; void FFT(Complex* a, int len, int DFT)//对a进行DFT或者逆DFT, 结果存在a当中,len必须是2的幂而且len大于多项式最高次数 { for(int i = 0; i < len; i++) A[rev(i, len)] = a[i]; for(int s = 1; (1 << s) <= len; s++) { int m = (1 << s); Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m)); for(int k = 0; k < len; k += m) { Complex w = Complex(1, 0); for(int j = 0; j < (m >> 1); j++) { Complex t = w*A[k + j + (m >> 1)]; Complex u = A[k + j]; A[k + j] = u + t; A[k + j + (m >> 1)] = u - t; w = w*wm; } } } if(DFT == -1) for(int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len; for(int i = 0; i < len; i++) a[i] = A[i]; return; } Complex a[maxn<<1]; int san[maxn]; int num[maxn<<1]; LL xishu[maxn<<1]; int main() { int t,n; cin>>t; while(t--){ scanf("%d",&n); int len=1; int s=0; for(int i=1;i<=n;i++){ scanf("%d",&san[i]); s=max(s,san[i]); } while(len<=s){ len<<=1; } len<<=1; sort(san+1,san+1+n); for(int i=0;i<len;i++){ a[i]=Complex(0,0); num[i]=0; } for(int i=1;i<=n;i++){ a[san[i]].real+=1; num[san[i]]++; } FFT(a,len,1); for(int i=0;i<len;i++){ a[i]=a[i]*a[i]; } FFT(a,len,-1); } }