题意：给出一堆数，问从这些数中取3个能组成三角形的概率？
sol：其实就是问从这些数里取3个组成三角形有多少种取法
脑洞大开的解法：用FFT
设一开始的数是1 3 3 4
作一个向量x，其中x[i]=边长为i的边的个数
那么就有x=[0 1 0 2 1 0 0 0 0]
令y=x，对x和y作DFT，得到dx和dy。令dn=dx*dy，再对dn作IDFT得到n
那么就得到n=[0 0 1 0 4 2 4 4 1 0 ]
其中n[i]=在x和y中各选一条边，使得两条边之和为i有几种方案
这时得到的n并不好，包含了各种重复的方案。还得减一下
最后得到的n=[0 0 0 0 2 1 1 2 0]
然后对n做奇怪的操作。。。
先求前缀和，得到S
对于a[i]. 我们假设a[i]是形成的三角形中最长的。这样就是在其余中选择两个和>a[i],而且长度不能大于a[i]的。（注意这里所谓的大于小于，不是说长度的大于小于，其实是排好序以后的，位置关系，这样就可以不用管长度相等的情况，排在a[i]前的就是小于的，后面的就是大于的）。 根据前面求得的结果。 长度和大于a[i]的取两个的取法是sum[len]-sum[a[i]]. 但是这里面有不符合的。 一个是包含了取一大一小的 cnt -= (long long)(n-1-i)*i; 一个是包含了取一个本身i,然后取其它的 cnt -= (n-1); 还有就是取两个都大于的了 cnt -= (long long)(n-1-i)*(n-i-2)/2;
Code：
1 #include <stdio.h> 2 #include <string.h> 3 #include <iostream> 4 #include <algorithm> 5 #include <math.h> 6 using namespace std; 7 const double PI = acos(-1.0); 8 #define LL long long 9 #define MAXL 524230 10 #define MAXN 100010 11 12 //复数结构体 13 struct Complex 14 { 15 double x,y;//实部和虚部 x+yi 16 Complex(double _x = 0,double _y = 0) 17 { 18 x = _x; 19 y = _y; 20 } 21 Complex operator -(const Complex &b)const 22 { 23 return Complex(x-b.x,y-b.y); 24 } 25 Complex operator +(const Complex &b)const 26 { 27 return Complex(x+b.x,y+b.y); 28 } 29 Complex operator *(const Complex &b)const 30 { 31 return Complex(x*b.x-y*b.y,x*b.y+y*b.x); 32 } 33 }; 34 /* 35 *进行FFT和IFFT前的反转变换。 36 *位置i和（i二进制反转后位置）互换 37 * len必须去2的幂 38 */ 39 void change(Complex y[],int len) 40 { 41 int i,j,k; 42 for(i = 1, j = len/2; i <len-1; i++) 43 { 44 if(i < j)swap(y[i],y[j]); 45 //交换互为小标反转的元素，i<j保证交换一次 46 //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的 47 k = len/2; 48 while(j >= k) 49 { 50 j -= k; 51 k /= 2; 52 } 53 if(j < k)j += k; 54 } 55 } 56 /* 57 *做FFT 58 * len必须为2^k形式， 59 * on==1时是DFT，on==-1时是IDFT 60 */ 61 //fft(x,len,1):对向量x做DFT（时域->频域），向量长度为1--len 62 //fft(x,len,-1):做IDFT（频域->时域） 63 void fft(Complex y[],int len,int on) 64 { 65 change(y,len); 66 for(int h = 2; h <= len; h <<= 1) 67 { 68 Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h)); 69 for(int j = 0; j < len; j+=h) 70 { 71 Complex w(1,0); 72 for(int k = j; k < j+h/2; k++) 73 { 74 Complex u = y[k]; 75 Complex t = w*y[k+h/2]; 76 y[k] = u+t; 77 y[k+h/2] = u-t; 78 w = w*wn; 79 } 80 } 81 } 82 if(on == -1) 83 for(int i = 0; i < len; i++) 84 y[i].x /= len; 85 } 86 87 Complex x[MAXL]; 88 LL nx[MAXL],tri[MAXN],sx[MAXL],NUM[MAXL]; 89 int TIMES,N,tm; 90 91 int main() 92 { 93 scanf("%d",&TIMES); 94 for (int Times=1;Times<=TIMES;Times++) 95 { 96 scanf("%d",&N); 97 98 memset(nx,0,sizeof(nx)); 99 memset(sx,0,sizeof(sx)); 100 memset(x,0,sizeof(x)); 101 memset(NUM,0,sizeof(NUM)); 102 103 for (int i=0;i<N;i++) 104 { 105 scanf("%d",&tri[i]); 106 NUM[tri[i]]++; 107 } 108 109 sort(tri,tri+N); 110 int mx=tri[N-1]; 111 int l1=mx+1; 112 int len=1; 113 while (len<2*l1) len<<=1; 114 for (int i=0;i<l1;i++) 115 x[i]=Complex(NUM[i],0); 116 for (int i=l1;i<len;i++) 117 x[i]=Complex(0,0); 118 119 fft(x,len,1); 120 121 for (int i=0;i<len;i++) 122 x[i]=x[i]*x[i]; 123 124 fft(x,len,-1); 125 for (int i=0;i<len;i++) 126 nx[i]=(LL)(x[i].x+0.5); 127 128 len=2*mx; 129 for (int i=0;i<N;i++) 130 nx[tri[i]*2]--; 131 for (int i=1;i<=len;i++) 132 nx[i]=nx[i]/2; 133 /* 134 for (int i=0;i<=len;i++) 135 cout<<nx[i]<<" "; 136 cout<<endl; 137 */ 138 sx[0]=0; 139 for (int i=1;i<=len;i++) 140 sx[i]=sx[i-1]+nx[i]; 141 /* 142 for (int i=0;i<=len;i++) 143 cout<<sx[i]<<" "; 144 cout<<endl; 145 */ 146 LL tot=0; 147 for (int i=0;i<N;i++) 148 { 149 tot+=sx[len]-sx[tri[i]]; 150 tot-=(LL)(N-i-1)*i; 151 tot-=(N-1); 152 tot-=(LL)(N-1-i)*(N-i-2)/2; 153 } 154 //cout<<tot<<endl; 155 LL kk=(LL)(N*(N-1)*(N-2)/6); 156 //cout<<tot<<" "<<kk<<endl; 157 printf("%.7lf ",(double)tot/kk); 158 } 159 160 return 0; 161 }
ref：http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html