HDU 4609 3-idiots(FFT学习)
2015-11-15 23:54
501 查看
附上kuangbin大神博客:讲解的非常清楚
分析:暴力枚举aka_k,那么问题变成了找满足ai+aj>ak(i,j,k互不相等)a_i+a_j>a_k(i,j,k互不相等)的(ai,aj)(a_i,a_j)的对数。FFT求出所有可能的和的对数,然后去掉坐标相同的就可以了。
PS:目前就写了两道FFT,基本上都是卷积的应用,FFT神奇的地就是能在nlognnlogn的时间内算出所有的X(k)=∑a[i]∗a[k−i]X(k)=\sum{a[i]*a[k-i]}
附上代码:
分析:暴力枚举aka_k,那么问题变成了找满足ai+aj>ak(i,j,k互不相等)a_i+a_j>a_k(i,j,k互不相等)的(ai,aj)(a_i,a_j)的对数。FFT求出所有可能的和的对数,然后去掉坐标相同的就可以了。
PS:目前就写了两道FFT,基本上都是卷积的应用,FFT神奇的地就是能在nlognnlogn的时间内算出所有的X(k)=∑a[i]∗a[k−i]X(k)=\sum{a[i]*a[k-i]}
附上代码:
#include <bits/stdc++.h> #define LL long long #define FOR(i,x,y) for(int i = x;i < y;++ i) #define IFOR(i,x,y) for(int i = x;i > y;-- i) using namespace std; //FFT copy from kuangbin const double pi = acos (-1.0); // Complex z = a + b * i struct Complex { double a, b; Complex(double _a=0.0,double _b=0.0):a(_a),b(_b){} Complex operator + (const Complex &rhs) const { return Complex(a + rhs.a , b + rhs.b); } Complex operator - (const Complex &rhs) const { return Complex(a - rhs.a , b - rhs.b); } Complex operator * (const Complex &rhs) const { return Complex(a * rhs.a - b * rhs.b , a * rhs.b + b * rhs.a); } }; //len = 2 ^ k void change (Complex y[] , int len) { for (int i = 1 , j = len / 2 ; i < len -1 ; i ++) { if (i < j) swap(y[i] , y[j]); int k = len / 2; while (j >= k) { j -= k; k /= 2; } if(j < k) j += k; } } // FFT // len = 2 ^ k // on = 1 DFT on = -1 IDFT 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].a /= len; } } } const int maxn = 100010; const int maxm = 100010; LL cnt[maxm<<2],sum[maxm<<2]; int a[maxn],n; Complex x[maxm<<2]; void work(){ memset(cnt,0,sizeof(cnt)); FOR(i,0,n) scanf("%d",&a[i]),++ cnt[a[i]]; sort(a,a+n); int mx_len = (a[n-1] + 1)<<1; int len = 1; while(len < mx_len) len <<= 1; FOR(i,0,len) x[i] = Complex(cnt[i],0.0); FFT(x,len,1); FOR(i,0,len) x[i] = x[i]*x[i]; FFT(x,len,-1); FOR(i,0,len) cnt[i] = (LL)(x[i].a+0.5); FOR(i,0,n) -- cnt[a[i]<<1]; FOR(i,0,len) cnt[i] >>= 1; sum[0] = cnt[0]; FOR(i,1,len) sum[i] = sum[i-1]+cnt[i]; LL ans = 0; FOR(i,0,n){ ans += sum[len-1] - sum[a[i]]; ans -= (LL)(n-i-1)*i; ans -= (LL)(n-1); ans -= (LL)(n-i-1)*(n-i-2)>>1; } LL tot = (LL)n*(n-1)*(n-2)/6; printf("%.7f\n",(double)ans/tot); } int main() { //freopen("test.in","r",stdin); int T; scanf("%d",&T); while(T--){ scanf("%d",&n); work(); } return 0; }
相关文章推荐
- 物联网与嵌入式系统概论-week2-Embedded Systems-Lesson3: Interacting with the Physical World
- 物联网与嵌入式系统概论-week2-Embedded Systems-Lesson2: Components of Embedded Systems
- 物联网与嵌入式系统概论-week2-Embedded Systems-Lesson1: Features and Constraints
- 物联网时代的35款开源工具
- 物联网版Win10 预览版10586 IoT Core版下载
- iOS - 音效 AudioToolbox.framework
- 【CTO讲堂】浅析工业级物联网项目的快速开发
- AUGTEK成为广域物联网国际联盟董事会唯一亚洲成员
- windows10 IOT +Azure会议概要总结
- windows10 IOT +Azure会议概要总结
- 初学嵌入式有感----------------推荐一个好老师及一个好的线上学习课程
- 树莓派从零安装物联网alljoyn环境
- 从DragonBoard看物联网技术
- IOS_AudioToolbox音效
- 浅谈工业级物联网项目架构设计及实施
- 物联网与嵌入式系统概论-week1-What Is the IoT- Lesson3: The Importance of the IoT
- 物联网与嵌入式系统概论-week1-What Is the IoT- Lesson2: Trends of the IoT
- 物联网与嵌入式系统概论-week1-What Is the IoT- Lesson1: Definition of the IoT
- IoTimerInLineHook
- nrf51822 GPIOTE