hdu_4609 3-idiots

题目:

  King OMeGa catched three men who had been streaking in the street. Looking as idiots though, the three men insisted that it was a kind of performance art, and begged the king to free them. Out of hatred to the real idiots, the king wanted to check if they were lying. The three men were sent to the king's forest, and each of them was asked to pick a branch one after another. If the three branches they bring back can form a triangle, their math ability would save them. Otherwise, they would be sent into jail.
  However, the three men were exactly idiots, and what they would do is only to pick the branches randomly. Certainly, they couldn't pick the same branch - but the one with the same length as another is available. Given the lengths of all branches in the forest, determine the probability that they would be saved.

思路:

  设 为被选出的三个数字(),那么可以列出下面三条式子:
  
  当 可以构成三角形时,这三条式子都成立。如果 不能构成三角形,因为 ,所以前两条式子成立,也就是说:
  
  所以只要求出 即可。利用生成函数, 对应的系数就是 等于 的方案数,减去相同木棍的方案。利用前缀和就能 算出成立的方案数了。总的复杂度

代码:

#include <bits/stdc++.h>
#define mset(a,b) memset(a,b,sizeof a)
#define _(...) (void)(__VA_ARGS__)
#define rep(i,a,b) for(int i(a),i##_END_(b);i<=i##_END_;i++)
#define drep(i,a,b) for(int i(a),i##_END_(b);i>=i##_END_;i--)
template<class T> inline bool tomax(T &a,T b){return a<b?a=b,1:0;}
template<class T> inline bool tomin(T &a,T b){return b<a?a=b,1:0;}
typedef long long ll;
const double PI=acos(-1);
const int N=262144+5;

struct Complex{
    double r,i;
    Complex operator + (const Complex _) const {
        return (Complex){r+_.r,i+_.i};
    }
    Complex operator - (const Complex _) const {
        return (Complex){r-_.r,i-_.i};
    }
    Complex operator * (const Complex _) const {
        return (Complex){r*_.r-i*_.i,i*_.r+_.i*r};
    }
};

int rev[N];
void fft(Complex a[],int n,bool inv){
    rep(i,0,n-1)
        if(rev[i]<i)
            std::swap(a[i],a[rev[i]]);
    for(int l=2,m=1;l<=n;l<<=1,m<<=1){
        Complex base=(Complex){cos(PI*2/l),(inv?-1:1)*sin(PI*2/l)};
        for(Complex *p=a;p!=a+n;p+=l){
            Complex w=(Complex){1,0};
            rep(i,0,m-1){
                Complex t=p[i+m]*w;
                p[i+m]=p[i]-t;
                p[i]=p[i]+t;
                w=w*base;
            }
        }
    }
    if(inv)
        rep(i,0,n-1)
            a[i].r/=n;
}

int log2(int x){
    int lim=0;
    while((1<<lim)<x)lim++;
    return lim;
}

Complex a[N];
ll cnt[N],sum[N],res[N];

int main(){
    int cas,n;
    scanf("%d",&cas);
    while(cas--){
        mset(cnt,0);mset(sum,0);mset(a,0);
        scanf("%d",&n);
        int mx=0,tmp;
        rep(i,1,n){
            scanf("%d",&tmp);
            a[tmp].r+=1;
            cnt[tmp]++;
            tomax(mx,tmp);
        }
        rep(i,1,mx*2)
            sum[i]=sum[i-1]+cnt[i];

        int lim=log2((mx+1)*2),d=1<<lim;
        rep(i,0,d-1)
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(lim-1));
        fft(a,d,false);
        rep(i,0,d-1)
            a[i]=a[i]*a[i];
        fft(a,d,true);

        rep(i,0,d-1){
            res[i]=a[i].r+0.5;
            if(i%2==0)
                res[i]-=cnt[i/2];
            res[i]/=2;
        }

        ll num=0;
        rep(i,0,mx*2)
            num+=res[i]*(sum[i-1]-2);
        ll all=(ll)n*(n-1)*(n-2)/6;
        ll ans=num-all*2;
        printf("%.7lf\n",(double)ans/all);
    }
    return 0;
}