BZOJ 3513: [MUTC2013]idiots

Description

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

Sol

FFT.

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

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

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

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

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

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
/**************************************************************
    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;
}

  

posted @   北北北北屿  阅读(148)  评论(0编辑  收藏  举报
编辑推荐:
· 深入理解 Mybatis 分库分表执行原理
· 如何打造一个高并发系统?
· .NET Core GC压缩(compact_phase)底层原理浅谈
· 现代计算机视觉入门之:什么是图片特征编码
· .NET 9 new features-C#13新的锁类型和语义
阅读排行:
· 手把手教你在本地部署DeepSeek R1,搭建web-ui ,建议收藏!
· Spring AI + Ollama 实现 deepseek-r1 的API服务和调用
· 《HelloGitHub》第 106 期
· 数据库服务器 SQL Server 版本升级公告
· C#/.NET/.NET Core技术前沿周刊 | 第 23 期(2025年1.20-1.26)
点击右上角即可分享
微信分享提示