BZOJ 3513: [MUTC2013]idiots

Description

求一个序列能够成的三角形个数.\(n \leqslant 10^5,a_i \leqslant 10^5\)

Sol

FFT.

我们可以用FFT求出任意两个形成的组合,不过要减去重复的.

我先算的是排列,最后除6变成组合.

然后考虑将第三条边加入,这时候只需要减去所有小于等于这条边的长度的个数*3即可.

因为这样正确的原因是我假设的一个条件,假设的是我们选的这条边是最长边,这样减掉的边都会比他短,也就考虑了所有方案.

因为算的是排列,并且其他两个的顺序已经确定,有三个位置可以选择,所以乘3.

Code

/**************************************************************
    Problem: 3513
    User: BeiYu
    Language: C++
    Result: Accepted
    Time:12652 ms
    Memory:36476 kb
****************************************************************/
 
#include <bits/stdc++.h>
using namespace std;
 
#define mpr make_pair
#define rr first
#define ii second
#define debug(a) cout<<#a<<"="<<a<<" "
typedef pair< double,double > Complex;
typedef long long LL;
const int N = 6e5+500;
const double Pi = M_PI;
 
Complex operator + (const Complex &a,const Complex &b) {
    return mpr(a.rr+b.rr,a.ii+b.ii);
}
Complex operator - (const Complex &a,const Complex &b) {
    return mpr(a.rr-b.rr,a.ii-b.ii);
}
Complex operator * (const Complex &a,const Complex &b) {
    return mpr(a.rr*b.rr-a.ii*b.ii,a.rr*b.ii+a.ii*b.rr);
}
inline int in(int x=0,char ch=getchar()) { while(ch>'9' || ch<'0') ch=getchar();
    while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();return x; }
 
 
int T,n,m;
LL ans;
Complex a[N],b[N],c[N];
LL w[N];
int l[N];
 
void init(int x) {
    for(n=1;n<=x;n<<=1);
    n<<=1;
}
void Rev(Complex a[]) {
    for(int i=0,j=0;i<n;i++) {
        if(i>j) swap(a[i],a[j]);
        for(int k=n>>1;(j^=k)<k;k>>=1);
    }
}
void DFT(Complex a[],int r=1) {
    Rev(a);
    for(int i=1;i<=n;i<<=1) {
        Complex wi=mpr(cos(2*Pi/i),r*sin(2*Pi/i));
        for(int k=0;k<n;k+=i) {
            Complex w=mpr(1.0,0.0);
            for(int j=k;j<k+i/2;j++) {
                Complex t1=a[j],t2=w*a[j+i/2];
                a[j]=t1+t2,a[j+i/2]=t1-t2;
                w=w*wi;
            }
        }
    }if(r==-1) for(int i=0;i<n;i++) a[i].rr/=n;
}
void FFT(Complex a[],Complex c[]) {
    DFT(a);
    for(int i=0;i<=n;i++) c[i]=a[i]*a[i];
    DFT(c,-1);
}
int main() {
    for(T=in();T--;) {
        memset(a,0,sizeof(a));
        m=n=in();
        int mx=0;
        for(int i=0;i<n;i++) {
            l[i]=in(),mx=max(mx,l[i]);
            a[l[i]].rr+=1;
        }
        if(m<3) { printf("%.7lf\n",0.0);continue; }
        init(mx);
        FFT(a,b);
        for(int i=0;i<n;i++) w[i]=b[i].rr+0.5;
         
        for(int i=0;i<m;i++) w[l[i]<<1]--;
//      for(int i=0;i<n;i++) cout<<w[i]<<" ";cout<<endl;
         
         
        for(int i=1;i<n;i++) w[i]+=w[i-1];
//      for(int i=0;i<n;i++) cout<<w[i]<<" ";cout<<endl;
         
        ans=(LL)m*(m-1)*(m-2);
        for(int i=0;i<m;i++) {
            ans-=w[l[i]]*3;
        }
         
        ans/=6;
        printf("%.7lf\n",(double)ans/((LL)m*(m-1)*(m-2)/6));    
    }
    return 0;
}