770
技術社區[雲棲]
HDU 4609 快速傅裏葉變換
題意：給出10^5個邊長，問任取三個組成三角形的概率。
這題沒遇到還不知道什麼是傅裏葉變換FFT，當時怎麼想時間複雜度都不能降到n*log(n)。實際是這樣的，把每個記錄長度個數的數組存在num數組中。
直接摘大神的吧 網址https://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html
學會了FFT。
這題就很容易了。
其實題目是給了n條線段。問隨機取三個，可以組成三角形的概率。
其實就是要求n條線段，選3條組成三角形的選法有多少種。
首先題目給了a數組，
如樣例一：
4
1 3 3 4
把這個數組轉化成num數組，num[i]表示長度為i的有num[i]條。
樣例一就是
num = {0 1 0 2 1}
代表長度0的有0根，長度為1的有1根，長度為2的有0根，長度為3的有兩根，長度為4的有1根。
使用FFT解決的問題就是num數組和num數組卷積。
num數組和num數組卷積的解決，其實就是從{1 3 3 4}取一個數，從{1 3 3 4}再取一個數，他們的和每個值各有多少個
例如{0 1 0 2 1}*{0 1 0 2 1} 卷積的結果應該是{0 0 1 0 4 2 4 4 1 }
長度為n的數組和長度為m的數組卷積，結果是長度為n+m-1的數組。
{0 1 0 2 1}*{0 1 0 2 1} 卷積的結果應該是{0 0 1 0 4 2 4 4 1 }。
這個結果的意義如下：
從{1 3 3 4}取一個數，從{1 3 3 4}再取一個數
取兩個數和為 2 的取法是一種：1+1
和為 4 的取法有四種：1+3， 1+3 ，3+1 ，3+1
和為 5 的取法有兩種：1+4 ，4+1；
和為 6的取法有四種：3+3,3+3,3+3,3+3,3+3
和為 7 的取法有四種： 3+4,3+4,4+3,4+3
和為 8 的取法有 一種：4+4
利用FFT可以快速求取循環卷積，具體求解過程不解釋了，就是DFT和FFT的基本理論了。
總之FFT就是快速求到了num和num卷積的結果。隻要長度滿足>=n+m+1.那麼就可以用循環卷積得到線性卷積了。
弄完FFT得到一個num數組，這個數組的含義在上麵解釋過了。
#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define N 400005
#define pi acos(-1.0) // PI值
struct complex
{
double r,i;
complex(double real=0.0,double image=0.0)
{
r=real;
i=image;
}
// 以下為三種虛數運算的定義
complex operator + (const complex o)
{
return complex(r+o.r,i+o.i);
}
complex operator - (const complex o)
{
return complex(r-o.r,i-o.i);
}
complex operator * (const complex o)
{
return complex(r*o.r-i*o.i,r*o.i+i*o.r);
}
} x1[N],x2[N];
void brc(complex *y,int l) // 二進製平攤反轉置換 O(logn)
{
register int i,j,k;
for(i=1,j=l/2; i<l-1; i++)
{
if(i<j) swap(y[i],y[j]); // 交換互為下標反轉的元素
// i<j保證隻交換一次
k=l/2;
while(j>=k) // 由最高位檢索，遇1變0，遇0變1，跳出
{
j-=k;
k/=2;
}
if(j<k) j+=k;
}
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
// 其中on==1時為DFT，on==-1為IDFT
{
register int h,i,j,k;
complex u,t;
brc(y,l); // 調用反轉置換
for(h=2; h<=l; h<<=1) // 控製層數
{
// 初始化單位複根
complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
for(j=0; j<l; j+=h) // 控製起始下標
{
complex w(1,0); // 初始化螺旋因子
for(k=j; k<j+h/2; k++) // 配對
{
u=y[k];
t=w*y[k+h/2];
y[k]=u+t;
y[k+h/2]=u-t;
w=w*wn; // 更新螺旋因子
} // 據說上麵的操作叫蝴蝶操作…
}
}
if(on==-1) for(i=0; i<l; i++) y[i].r/=l; // IDFT
}
long long num[N];
long long sum[N];
int a[N>>1];
int main()
{
int t,n,l,maxn;
scanf("%d",&t);
while(t--)
{
memset(num,0,sizeof(num));
maxn=0;
scanf("%d",&n);
for(int i=0; i<n; i++)
{
scanf("%d",&a[i]);
num[a[i]]++;
}
sort(a,a+n);
maxn=a[n-1];
l=1;
while(l<2*(maxn+1)) l<<=1;
for(int i=0; i<=maxn; i++)
x1[i].r=num[i],x1[i].i=0;
for(int i=maxn+1; i<=l; i++)
x1[i].r=x1[i].i=0;
fft(x1,l,1);
for(int i=0; i<=l; i++) x1[i]=x1[i]*x1[i];
fft(x1,l,-1);
for(int i=0; i<=2*maxn; i++) num[i]=(long long)(x1[i].r+0.5);
for(int i=0; i<n; i++) num[a[i]+a[i]]--;
for(int i=1; i<=2*maxn; i++) num[i]/=2;
sum[0]=0;
for(int i=1; i<=2*maxn; i++)
sum[i]=sum[i-1]+num[i];
long long cnt=0;
for(int i=0; i<n; i++) //枚舉每條邊並把這條邊當成最大邊
{
cnt+=sum[2*maxn]-sum[a[i]];
cnt-=n-1; //減去包括這條邊的
cnt-=(long long)i*(n-i-1);//減去有一條邊比這條邊大的情況
cnt-=(long long)(n-i-1)*(n-i-2)/2;//減去有兩條都比該條邊大的情況
}
long long tol=(long long)n*(n-1)*(n-2)/6;
printf("%.7f\n",double((double)cnt/(double)tol));
}
return 0;
}
最後更新：2017-04-03 16:48:43