hdu 4609 3-idiots <FFT>

链接: http://acm.hdu.edu.cn/showproblem.php?pid=4609

题意: 给定 N 个正整数, 表示 N 条线段的长度, 问任取 3 条, 可以构成三角形的概率为多少~

数据范围: N<=10^5 ~~

思路:设三边分别为 x, y, z (x<=y<=z) 枚举 z ,统计 x+y 大于 z 的数目 .

比赛时能想到的只有 O(n^2) 的算法,无力 AC~

赛后才知道有种东西叫 FFT ~

以下为官方解题报告:

/*

　　记录 A_i 为长度为 i 的树枝的数量，并让 A 对它本身做 FFT，得到任意选两个树枝能得到的各个和的数量。枚举第三边，

计算出所有两边之和大于第三条边的方案数，并把前两条边包含最长边的情况减掉就是答案。

*/

设长度为 a1,a2, ....an, 统计每种长度个数为 cnt[ ai ] ; 可表示多项式      

多项式自乘后， 其指数即为 ai + aj 的值， 系数即为 方案数， 减去不合法的即为所求~

#include <cstdio>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;
typedef __int64 LL;
const double pi=acos(-1);
const int MAXN=3e5;
const double eps=1e-6;
struct C
{
    double r, i;
    C (){}
    C(double _r, double _i):r(_r),i(_i){}
    inline C operator +(const C a)const{
    return C(r+a.r, i+a.i);    
    }
    inline C operator - (const C a)const {
        return C(r-a.r, i-a.i);    
    }
    inline C operator * (const C a)const{
        return C(r*a.r-i*a.i, r*a.i+i*a.r);    
    }
}a[MAXN], b[MAXN];
int num[MAXN], cnt[MAXN];
LL res[MAXN], sum[MAXN];

void brc(C *y,int l) // 二进制平摊反转置换 O(logn)
{
    register int i,j,k;
    for(i=1,j=l>>1;i<l-1;i++){
        if(i<j)    swap(y[i],y[j]); // 交换互为下标反转的元素
                                // i<j保证只交换一次
        k=l>>1;
        while(j>=k) {// 由最高位检索，遇1变0，遇0变1，跳出
            j-=k;
            k>>=1;
        }
        if(j<k)    j+=k;
    }
}
void FFT(C *y,int l,int on) // FFT O(nlogn)
                             // 其中on==1时为DFT，on==-1为IDFT
{
    register int h,i,j,k;
    C u,t; 
    brc(y,l); // 调用反转置换
    for(h=2;h<=l;h<<=1) // 控制层数
    {
        // 初始化单位复根
        C wn(cos(on*2*pi/h),sin(on*2*pi/h));
        for(j=0;j<l;j+=h) // 控制起始下标
        {
            C 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
}

int n,  L , Max, T;
int main() {
    scanf("%d", &T);
    while (T--) {
        Max = 0;
        memset(cnt, 0, sizeof(cnt));
        scanf("%d", &n);
        for (int i = 0; i < n; ++i) {
            scanf("%d", num + i);
            cnt[num[i]]++;
            Max = max(Max, num[i]);
        }
        ++Max;L = 1;
        while (L < Max <<1) {
            L <<= 1;
        }
        for (int i = 0; i < Max; ++i) {
            a[i] = C(cnt[i], 0);
        }

        for (int i = Max; i < L; ++i) {
            a[i] = C(0, 0);
        }
        FFT(a, L, 1);
        for (int i = 0; i < L; ++i)
            a[i] = a[i] * a[i]; // 多项式自乘
        FFT(a, L, -1);

        for (int i = 0; i < L; ++i) {
            res[i] = (LL) (a[i].r + 0.5);
        }

        for (int i = 0; i <= Max; ++i) 
            res[i << 1] -= cnt[i];// 自己和自己乘， 即 枚举的是 x+x 的 
    
        for (int i = 0; i < L; ++i)
            res[i] >>= 1; //  x+y 与  y+x 算一次 
        for (int i = 1; i < L; ++i) {  //前缀和 
            sum[i] = sum[i - 1] + res[i];     
        }
        double tot = 0, den =  1.0*n * (n - 1) * (n - 2)/6;
        
        for (int i = 0; i < n; ++i) {   
            tot += sum[num[i]] / den;//  x+y<=z 的个数 
        
        }        
        double ans = 1 - tot ;
        printf("%.7f\n", ans);
    }

    return 0;
}