首页 > 代码库 > 快速傅里叶变换FFT——kuangbin模板
快速傅里叶变换FFT——kuangbin模板
/* #include <bits/stdc++.h> using namespace std; const double PI = acos(-1.0); struct Complex { double r,i; Complex(double _r = 0,double _i = 0) { r = _r; i = _i; } Complex operator +(const Complex &b) { return Complex(r+b.r,i+b.i); } Complex operator -(const Complex&b) { return Complex(r-b.r,i-b.i); } Complex operator *(const Complex &b) { return Complex(r*b.r-i*b.i,r*b.i+i*b.r); } }; void change(Complex y[],int len) { int i,j,k; for(i = 1, j = len/2;i < len-1;i++) { if(i < j)swap(y[i],y[j]); k = len/2; while( j >= k) { j -= k; k /= 2; } if(j < k)j += k; } } void fft(Complex y[],int len,int on) { change(y,len); for(int h = 2;h <= len;h <<= 1) { Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h)); for(int j = 0;j < len;j += h) { Complex w(1,0); for(int k = j;k < j+h/2;k++) { Complex u = y[k]; Complex t = w*y[k+h/2]; y[k] = u+t; y[k+h/2] = u-t; w = w*wn; } } } if(on == -1) for(int i = 0;i < len;i++) y[i].r /= len; } const int MAXN = 400040; Complex x1[MAXN]; int a[MAXN/4]; long long num[MAXN];//100000*100000会超int long long sum[MAXN]; int main() { int T; int n; scanf("%d",&T); while(T--) { scanf("%d",&n); memset(num,0,sizeof(num)); for(int i = 0;i < n;i++) { scanf("%d",&a[i]); num[a[i]]++; } sort(a,a+n); int len1 = a[n-1]+1; int len = 1; while( len < 2*len1 )len <<= 1; for(int i = 0;i < len1;i++) x1[i] = Complex(num[i],0); for(int i = len1;i < len;i++) x1[i] = Complex(0,0); fft(x1,len,1); for(int i = 0;i < len;i++) x1[i] = x1[i]*x1[i]; fft(x1,len,-1); for(int i = 0;i < len;i++) num[i] = (long long)(x1[i].r+0.5); len = 2*a[n-1]; //减掉取两个相同的组合 for(int i = 0;i < n;i++) num[a[i]+a[i]]--; //选择的无序,除以2 for(int i = 1;i <= len;i++) { num[i]/=2; } sum[0] = 0; for(int i = 1;i <= len;i++) sum[i] = sum[i-1]+num[i]; long long cnt = 0; for(int i = 0;i < n;i++) { cnt += sum[len]-sum[a[i]]; //减掉一个取大,一个取小的 cnt -= (long long)(n-1-i)*i; //减掉一个取本身,另外一个取其它 cnt -= (n-1); //减掉大于它的取两个的组合 cnt -= (long long)(n-1-i)*(n-i-2)/2; } //总数 long long tot = (long long)n*(n-1)*(n-2)/6; printf("%.7f\n",(double)cnt/tot); } return 0; }*/ #include <bits/stdc++.h> using namespace std; const double PI = acos(-1.0); const int MAXN = 400044; struct Complex{ double x,y; Complex(double _x = 0.0,double _y =0.0){ x = _x;y =_y; } Complex operator - (const Complex &b)const{ return Complex(x-b.x,y-b.y); } Complex operator + (const Complex &b)const{ return Complex(x+b.x,y+b.y); } Complex operator * (const Complex &b)const{ return Complex(x*b.x - y*b.y,x*b.y + y*b.x); } }; void change (Complex y[],int len){ int i,j,k; for (i = 1,j = len/2;i<len -1;i++){ if (i<j) swap(y[i],y[j]); k = len/2; while (j >= k){ j-=k; k/=2; } if (j<k) j+=k; } } void fft(Complex y[],int len, int on) { change(y,len); for (int h = 2; h<=len;h<<=1){ Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h)); for (int j =0 ; j< len; j+=h){ Complex w(1,0); for (int k = j;k < j+h/2;k++){ Complex u = y[k]; Complex t = w*y[k+h/2]; y[k] = u+t; y[k+h/2] = u-t; w = w*wn; } } } if (on == -1) for (int i = 0; i<len; i++) y[i].x /= len; } Complex x1[MAXN]; int a[MAXN/4]; long long num[MAXN]; long long sum[MAXN]; int main() { int T; cin >> T; int n; while (T--){ scanf("%d",&n); memset(num,0,sizeof(num)); for (int i = 0; i< n; i++){ scanf("%d",&a[i]); num[a[i]]++; } sort(a,a+n); int len1 = a[n-1]+1; int len = 1; while ( len < 2*len1 ) len <<= 1; for (int i = 0; i<len1 ; i++) x1[i] = Complex(num[i],0); for (int i = len1 ; i < len ;i++) x1[i] = Complex(0,0); fft(x1,len,1); for (int i = 0; i< len; i++) x1[i] = x1[i] * x1[i]; fft(x1,len,-1); for (int i = 0 ; i<len ; i++) num[i] = (long long) (x1[i].x + 0.5); len = 2*a[n-1]; for (int i = 0 ; i<n; i++) num[a[i]+a[i]]-- ; for (int i = 1; i<=len;i++) num[i]/=2; sum[0] = 0; for (int i = 1; i<=len ; i++) sum[i] = sum[i-1] + num[i]; long long cnt = 0; for (int i = 0; i <n; i++){ cnt += sum[len] - sum[a[i]]; cnt -= (long long)(n-1-i)*i; cnt -= (n-1); cnt -= (long long)(n-1-i)*(n-i-2)/2; } long long tot = (long long) n * (n-1)*(n-2)/6; printf("%.7f\n",(double)(cnt*1.0/tot)); } return 0; }
快速傅里叶变换FFT——kuangbin模板
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。