HDU 4609 3-idiots FFT入门

最新推荐文章于 2020-10-08 21:17:34 发布

npugenji

最新推荐文章于 2020-10-08 21:17:34 发布

阅读量218

点赞数

CC 4.0 BY-SA版权

分类专栏： FFT

本文链接：https://blog.youkuaiyun.com/z631681297/article/details/78128473

FFT 专栏收录该内容

4 篇文章

订阅专栏

本文介绍了一种使用快速傅立叶变换（FFT）计算由一组数构成三角形概率的方法。通过计算两数之和并考虑重复计数的问题，最终得出正确概率。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

给你一堆数，问你组成三角形的概率是多少。

首先用FFT算出可以用两个数组成的长度。

把数当做权值，出现一个数就在那个权值位置加1。
然后就是用FFT计算多项式乘法。

然而这样会有重复
1. 一个数自己和自己组合是不存在的，所以要删一次
2. 比如有两个数1 2 ，取(1,2)和取(2,1)是一样的，这样等于所有情况都算了两遍，除以2即可。

然后算能组成三角形的情况个数。
把FFT的结果先做一个前缀和。
为了避免重复，把原来的数组排序，默认我们取的那根长度是三角形里最长的。
1. 那么组合出的长度比自己长的情况有sum[ maxlen*2 ] - sum[ a[i] ] 种。
2. 取的木棒不能再和其他的组合了减去n-1种情况。
3. 不能用比自己长的木棒，减去 $C_{i-1}^{2}$
4. 不能一个比自己长，一个比自己短，减去 $(i-1)*(n-i)$

最后除 $C_{n}^{3}$ 就是答案了

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include<cmath>

using namespace std;
typedef long long LL;
const int INF = 0x3f3f3f3f;
const LL INFL = 0x3f3f3f3f3f3f3f3fLL;
const double eps = 1e-8;
const double pi = acos(-1.0);
const int MAXN = 131072*4+5;


struct cp
{
    double x,y;
    cp() {}
    cp(double x,double y):x(x),y(y) {}
    inline double real() { return x; }
    inline cp operator * (const cp& r) const { return cp(x*r.x -y*r.y,x*r.y+y*r.x);  }
    inline cp operator - (const cp& r) const { return cp(x-r.x,y-r.y); }
    inline cp operator + (const cp& r) const { return cp(x+r.x,y+r.y); }
};

cp a[MAXN],b[MAXN];
LL r[MAXN],res[MAXN];
LL ax[MAXN],l[MAXN];
//ax bx res 应该是同一类型,要开就都开LL

void fft_init(int nm,int k)
{
    for (int i=0;i<nm;i++) r[i] = (r[i>>1] >>1) | ((i & 1) << (k-1));
}

void fft(cp ax[],int nm,int op)
{
    for (int i=0;i<nm;i++) if (i<r[i]) swap(ax[i],ax[r[i]]);
    for (int h=2,m=1;h<=nm;h<<=1,m<<=1)
    {
        cp wn = cp(cos(op*2*pi/h),sin(op*2*pi/h));
        for (int i=0;i<nm;i+=h)
        {
            cp w(1,0);
            for (int j=i;j<i+m;++j,w=w*wn)
            {
                cp t=w*ax[j+m];
                ax[j+m] = ax[j] -t;
                ax[j] = ax[j] +t;
            }
        }
    }
    if (op == -1) for (int i=0;i<nm;i++) ax[i].x /= nm;
}

void trans(LL ax[],LL bx[],int n,int m)
{
    int nm=1,k=0;
    while (nm < 2*n || nm<2*m) nm<<=1,++k;

    for (int i=0;i<n;i++) a[i] = cp(ax[i],0);
    for (int i=0;i<m;i++) b[i] = cp(bx[i],0);
    for (int i=n;i<nm;i++) a[i] = cp(0,0);
    for (int i=m;i<nm;i++) b[i] = cp(0,0);

    fft_init(nm,k);
    fft(a,nm,1);fft(b,nm,1);
    for (int i=0;i<nm;i++) a[i] = a[i] * b[i];
    fft(a,nm,-1);
    nm = n+m-1;
    for (int i=0;i<nm;i++)
    {
        res[i] = (LL)(a[i].real() + 0.5);
    }
}

int n,m;


bool cmp(int a,int b)
{
    return a>b;
}

int main()
{
    int t;
    scanf("%d",&t);
    while (t--)
    {
        memset(res,0,sizeof res);
        memset(ax,0,sizeof ax);
        LL anum = 0;
        LL ans=0;
        scanf("%d",&n);
        LL tot = 1ll*n*(n-1)*(n-2)/6;
        for (int i=1;i<=n;i++)
        {
            scanf("%d",&l[i]);
            ax[l[i]] ++;
            anum = max(anum,l[i]);
        }
        trans(ax,ax,anum+1,anum+1);
        sort(l+1,l+1+n,cmp);
        for (int i=1;i<=n;i++) res[ 2*l[i] ] --;
        for (int i=1;i<=anum*2;i++) res[i] /= 2 ;//选1 2 和选 2 1 是一样的
        for (int i=0;i<=anum*2;i++) res[i] += res[i-1];
        for (int i=1;i<=n;i++)
        {
            ans += res[anum*2] - res[l[i]];
            ans -= 1ll*(i-1)*(i-2)/2;//两个大
            ans -= 1ll*(i-1)*(n-i);//一大一小
            ans -= 1ll*(n-1);//自己和别人
        }
        printf("%.7f\n",1.0*ans/tot);
    }
    return 0;
}