HDU 4609 3-idiots FFT入门

来源:互联网 发布:山西软件开发公司 编辑:程序博客网 时间:2024/06/01 09:45

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

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

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

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

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

最后除C3n就是答案了

#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 应该是同一类型,要开就都开LLvoid 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;}
原创粉丝点击