HDU 4609-526互联

description

给定一个长度为 \(n\) 的序列 \(A\)，元素值域大小为 \(10^5\)。求从中任选三个不同位置的元素，以它们的值为三边能够成三角形的概率。

solution

设有 \(cnt\) 种选三个不同的元素构成三角形的方案，则答案显然为 \(\dfrac{6cnt}{n(n-1)(n-2)}\)。

\(a,b,c(a\leq b\leq c)\) 能作为三角形三边的充要条件是 \(a+b>c\)。但是由于要求 \(a\leq b\leq c\)，枚举 \(a+b\)，我们不好直接统计出有多少个 \(c\) 满足偏序关系且 \(a+b>c\)。

但是可以反着来，统计有多少对 \((a,b,c)\) 构不成三角形，也就是 \(a+b\leq c\)。这就很好统计了。预处理 \(cc_p\) 表示大于等于 \(p\) 的元素的数量，然后构造多项式统计 \(a+b=p\) 的数量 \(x_p\)，枚举 \(p\)，把 \(cc_px_p\) 贡献到答案。输出时记得把答案变成合法的概率。

code

#include<bits/stdc++.h>

using namespace std;

const int N=(1<<22)+10;

typedef complex<double> E;
const double pi=acos(-1);

int lstn,rev[N];
void prework(int n){
  if(lstn==n) return ;
  for(int i=1; i<n; i++) rev[i]=(rev[i^(i&-i)]|(1<<(__lg(n)-__lg(i&-i)-1)));
  return lstn=n,void();
}

void fft(E *p,int n,int op){
  prework(n);
  for(int i=1; i<n; i++) if(i<rev[i]) swap(p[i],p[rev[i]]);

  for(int i=2; i<=n; i<<=1){
    int len=i>>1;
    E w=cos(2*pi/i)+1i*sin(2*pi/i);
    if(op==-1) w-=2i*sin(2*pi/i);
    for(int pos=0; pos<n; pos+=i){
      E now=1;
      for(int j=pos; j<pos+len; j++,now=now*w){
        auto u=p[j],v=p[j+len]*now;
        p[j]=u+v,p[j+len]=u-v;
      }
    }
  }

  if(op==-1){
    for(int i=0; i<n; i++) p[i]/=n;
  }
}

E a[N];
long long n,cnt1,s[N],maxa,cnt2,cc[N];

void solve(){
  cin>>n;
  cnt2=0;
  cnt1=n*(n-1)*(n-2)/6;
  maxa=0;
  for(int i=1; i<=n; i++){
    int x;
    scanf("%d",&x);
    s[x]++; cc[x]++;
    a[x]+=1.0;
    maxa=max(maxa,x*1ll);
  }

  for(int i=maxa-1; ~i; i--) s[i]+=s[i+1];

  int t=1;
  for(;t<=maxa*2;t<<=1) ;

  fft(a,t,1);
  for(int i=0; i<t; i++) a[i]=a[i]*a[i];
  fft(a,t,-1);

  for(int i=0; i<t; i++){
    long long x=(long long)(a[i].real()+0.5);
    if(i%2==0){
       x-=(cc[i/2]);
    }
    x/=2;
    cnt2+=s[i]*x;
  }

  printf("%.7lf\n",((double)(cnt1-cnt2))/cnt1);

  memset(s,0,sizeof(long long)*(maxa+10));
  memset(cc,0,sizeof(long long)*(maxa+10));
  memset(a,0,sizeof(E)*(t+1));
}

int main(){

  int T;
  cin>>T;
  while(T--) solve();

  return 0;
}