bzoj3513 [MUTC2013]idiots（生成函数+FFT）

题目链接

分析：
我们受到bzoj3771的启发
可以用生成函数计算边之和为n的方案数
设计生成函数$A(x)$，如果有长度为n的木棍，$x^n$的系数就+1

如何判断三边是否成三角形呢？
著名的三角不等式：两边之和大于第三边

那就可以衍生出一种方法：
先算出选出两条边，长度和为n的方案数
那么第三条边一定小于n，累计长度小于n的边数，乘上$x^n$的系数即可

但是这样存在一个小问题：我们选择的第三条边，有可能和长度和为n的那两条边重复
我们可以用容斥原理处理一下吗？
很遗憾，并不能
因为我们并不清楚“长度和为n”的组合情况是什么，不好直接容斥

怎么破？

想起我们的原则：正难则反

既然我们直接计算答案比较困难，就可以考虑计算答案的补集

下面我们就计算两边之和小于等于第三边的方案数：

两边之和为X的方案数：
$A(x)*A(x)=A^2(x)$
但是其中有同一个木棍选了两次的情况，我们需要减掉
设计$B(x)$，表示强制一个木棍选了两次的情况：如果有长度为n的木棍，$x^{2n}$的系数就+1
不重复的选择两条木棍的方案数：$A^2(x)-B(x)$
但是这样得到的是两条木棍的排列数，还要除以2：$C(x)={A^2(x)-B(x) \over 2}$

我们将$C(x)$求一个前缀和，表示两边之和小于等于X的方案数
$A(x)*C(x)$即为所求
因为$C$中是两边之和小于等于X的方案数，那么我们选择的第三边一定大于那两条边，不用去重了

因为计算的是概率，最后不要忘了除以总方案数

$ans={1-{A(x)*C(x) \over C(n,3)}}$

tip

注意两边之和的范围

为了方便计算，C和B直接用实数计算即可

#include<bits/stdc++.h>
#define ll long long

using namespace std;

const double pi=acos(-1.0);
const int N=300010;
struct node{
    double x,y;
    node (double xx=0,double yy=0) {
        x=xx,y=yy;
    }
};
node A[N],o[N],_o[N];
int n,fn,B[N];
ll C[N];

node operator +(const node &a,const node &b) {return node(a.x+b.x,a.y+b.y);}
node operator -(const node &a,const node &b) {return node(a.x-b.x,a.y-b.y);}
node operator *(const node &a,const node &b) {return node(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

void init(int n) {
    for (int i=0;i<=n;i++) {
        o[i]=node(cos(2.0*i*pi/n),sin(2.0*i*pi/n));
        _o[i]=node(cos(2.0*i*pi/n),-sin(2.0*i*pi/n));
    }
}

void FFT(int n,node *a,node *w) {
    int i,j=0,k;
    for (i=0;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>1;(j^=l)<l;l>>=1);
    }
    for (i=2;i<=n;i<<=1) {    //i:合并的区间长度 
        int m=i>>1;
        for (j=0;j<n;j+=i) 
            for (k=0;k<m;k++) {
                node z=a[j+m+k]*w[n/i*k];   //
                a[j+m+k]=a[j+k]-z;
                a[j+k]=a[j+k]+z;
            }
    }
}

int main()
{
    int T;
    scanf("%d",&T);
    while (T--) {
        for (int i=0;i<=fn;i++) {
            A[i].x=0; A[i].y=0;
            B[i]=0;
        }

        scanf("%d",&n);
        int mx=0;
        for (int i=0;i<n;i++) {
            int x;
            scanf("%d",&x);
            A[x].x++; B[x]++;
            mx=max(mx,x);
        }

        fn=1;
        while (fn<=mx+mx) fn<<=1;
        init(fn);
        FFT(fn,A,o);
        for (int i=0;i<=fn;i++) 
            A[i]=A[i]*A[i];
        FFT(fn,A,_o);
        for (int i=0;i<=fn;i++)   //A*A-B 
        {
            C[i]=(ll)(A[i].x/fn+0.5);
            if (i%2==0) C[i]-=B[i>>1];   //长度偶数才会出现重复选择的情况 
        }
        for (int i=1;i<=fn;i++)   //前缀和 
            C[i]=C[i]+C[i-1];

        double ans=0;
        double c=(double)n*(n-1)*(n-2)/6.0;
        for (int i=1;i<=mx;i++)   //注意两边之和的范围 
            ans+=(double)C[i]*B[i];
        ans/=2.0;                 //组合 
        ans=1.0-ans/c;
        printf("%0.7lf\n",ans);
    }
}