题目链接：https://vjudge.net/problem/HDU-4609
题目大意：
给出 \(n\) 条线段长度，问任取 \(3\) 根，组成三角形的概率。
知识点： FFT、生成函数
解题思路：
\(aX^b\) 代表长度为 \(b\) 的线段有 \(a\) 条，统计各种长度的线段的条数，列出多项式：\(\sum a_{i} X^{i}\)，其中 \(a_i\) 代表长度为 \(i\) 的线段的条数。
则由 \(\sum a_{i} X^{i} * \sum a_{i} X^{i}\) 得出另一条多项式：\(\sum b_{i} Y^{i}\)，其中 \(b_{i}\) 代表从 \(n\) 条线段中有放回地取两条线段，长度和为 \(i\) 的方案数。显然，这些方案中有不合法的方案：一种情况是同一条线段取两次；另一种情况是两条不同的线段由于取的次序不同被视为两种方案，在此处应该进行第一次去重。最后计算前缀和：\(sum(i)\) 表示长度大于等于 \(i\) 的方案数。
接下来求能组成三角形的方案数。对线段从短到长进行排序（注：下标从 \(0\) 开始）。每一条线段对于能组成三角形的方案数的贡献即为找出两条下标比它小并且加起来总长度比它大的方案数。先加上 \(sum(i+1)\)，但这里面也有不合法的：一种是这条线段本身再加上另一条线段（减去\((N-1)\)）；一种是两条下标比它大的线段（减去\(C_{N-1-i}^2\)）;最后一种是选择的两条线段中一条下标比它小，一条比它大（减去 \(i(N-1-i)\) ).
最终答案即为能组成三角形的方案数除以总的取线方案\((C_N^3)\).
写到这里，忽然想起中二帽曾经说过的一句话：“多项式乘法本质上就是一种卷积。”
AC代码：
1 #include <bits/stdc++.h> 2 3 using namespace std; 4 typedef long long ll; 5 const double PI = acos(-1.0); 6 const int maxn = 1e5+5; 7 ll num[maxn]; 8 struct Complex{ 9 double x,y; 10 Complex(double _x = 0.0, double _y = 0.0){ 11 x = _x; 12 y = _y; 13 } 14 Complex operator -(const Complex &b)const{ 15 return Complex(x-b.x,y-b.y); 16 } 17 Complex operator +(const Complex &b)const{ 18 return Complex(x+b.x,y+b.y); 19 } 20 Complex operator *(const Complex &b)const{ 21 return Complex(x*b.x-y*b.y,x*b.y+y*b.x); 22 } 23 }; 24 void change(Complex y[],int len){ 25 int i,j,k; 26 for(i=1,j=len/2;i<len-1;i++){ 27 if(i<j) swap(y[i],y[j]); 28 k=len/2; 29 while(j>=k){ 30 j-=k; 31 k/=2; 32 } 33 if(j<k) j+=k; 34 } 35 } 36 void fft(Complex y[],int len,int on){ 37 change(y,len); 38 for(int h=2;h<=len;h<<=1){ 39 Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h)); 40 for(int j=0;j<len;j+=h){ 41 Complex w(1,0); 42 for(int k=j;k<j+h/2;k++){ 43 Complex u=y[k]; 44 Complex t=w*y[k+h/2]; 45 y[k]=u+t; 46 y[k+h/2]=u-t; 47 w=w*wn; 48 } 49 } 50 } 51 if(on==-1) 52 for(int i=0;i<len;i++) 53 y[i].x/=len; 54 } 55 Complex x[maxn<<2]; 56 ll have[maxn<<2],sum[maxn<<2]; 57 int a[maxn]; 58 int main(){ 59 int T,N; 60 scanf("%d",&T); 61 while(T--){ 62 memset(num,0,sizeof(num)); 63 scanf("%d",&N); 64 for(int i=0;i<N;i++){ 65 scanf("%d",&a[i]); 66 num[a[i]]++; 67 } 68 sort(a,a+N); 69 int limt=a[N-1]+1; 70 int len=1; 71 while(len<limt*2) len<<=1; 72 for(int i=0;i<=limt;i++) x[i]=Complex(num[i],0); 73 for(int i=limt+1;i<len;i++) x[i]=Complex(0,0); 74 fft(x,len,1); 75 for(int i=0;i<len;i++) 76 x[i]=x[i]*x[i]; 77 fft(x,len,-1); 78 for(int i=0;i<len;i++) 79 have[i]=(ll)(x[i].x+0.5); 80 len=2*a[N-1]; 81 for(int i=0;i<N;i++) 82 have[a[i]*2]--; 83 for(int i=1;i<=len;i++) have[i]/=2; 84 85 sum[len+1]=0; 86 for(int i=len;i>=0;i--) 87 sum[i]=sum[i+1]+have[i]; 88 sort(a,a+N); 89 ll ans=0; 90 for(int i=0;i<N;i++){ 91 ans+=sum[a[i]+1]; 92 ans-=(N-1); 93 ans-=(ll)i*(N-1-i); 94 ans-=(ll)(N-1-i)*(N-2-i)/2; 95 } 96 double tot=(double)N*(N-1)*(N-2)/6.0; 97 printf("%.7lf\n",(double)ans/tot); 98 } 99 return 0; 100 }
