HDU 4602 FFT 终于大致会用FFT了。。

73 篇文章 21 订阅
56 篇文章 0 订阅


还是去看他的题解吧。。好详细的说。。kuangbin大牛啊  


http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210565.html


直接贴代码。。

/*
 *@author ipqhjjybj
 *@date 2013/08/12
 */
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
const int N=1<<18;
#define ll long long
const double pi=acos(-1.0);
using namespace std;

struct Complex{
    double r,i;
    Complex(double real=0.0,double image=0.0){
        r=real;i=image;
    }
    Complex operator + (const Complex &o){
        return Complex(r+o.r,i+o.i);
    }
    Complex operator - (const Complex &o){
        return Complex(r-o.r,i-o.i);
    }
    Complex operator * (const Complex &o){
        return Complex(r*o.r-i*o.i,r*o.i+i*o.r);
    }
}x1[N],x2[N];
int a[N];
ll sum[N],p[N];
int l;
void brc(Complex *y,int l){
    register int i,j,k;
    for(i=1,j=(l>>1);i<l-1;i++){
        if(i<j) swap(y[i],y[j]);
        k=l>>1;
        while(j>=k){
            j-=k;
            k/=2;
        }
        if(j<k) j+=k;
    }
}
void fft(Complex *y,int l,double on){
    register int h,i,j,k;
    Complex u,t;
    brc(y,l);
    for(h=2;h<=l;h<<=1){
        Complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
        for(j=0;j<l;j+=h){
            Complex w(1,0);
            for(k=j;k<j+(h>>1);k++){
                u=y[k];
                t=w*y[k+(h>>1)];
                y[k]=u+t;
                y[k+(h>>1)]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1) for(i=0;i<l;i++)y[i].r/=l;
}
int main(void){
    //freopen("4609.in","r",stdin);
    register int i;
    int t,n;
    scanf("%d",&t);
    while(t--){
        l=1<<19;
        scanf("%d",&n);
        for(i=0;i<n;i++)
            scanf("%d",a+i);
        sort(a,a+n);
        while((l>>2)>a[n-1]) l>>=1;
        for(i=0;i<l;i++)
            sum[i]=p[i]=0;
        for(i=0;i<n;i++)
            p[a[i]]++;
        for(i=0;i<l;i++){
            x2[i]=x1[i]=Complex(1.0*p[i],0.0);
        }
        fft(x1,l,1);
        fft(x2,l,1);
        for(i=0;i<l;i++)x1[i]=x1[i]*x2[i];
        fft(x1,l,-1);
        for(i=0;i<l;i++)sum[i]=(ll)(x1[i].r+0.5);

        //for(i=0;i<l;i++){
        //    printf("sum[%d]=%lld\n",i,sum[i]);
        //}

        for(i=0;i<n;i++)sum[a[i]+a[i]]--;
        for(i=1;i<l;i++)sum[i]>>=1;
        ll tot=(ll)n*(n-1)*(n-2)/6;
        ll ans=tot;
        int k=0;
        for(i=0;i<l;i++)
            if(sum[i]){
                while(k<n&&a[k]<i)k++;
                if(k==n)break;
                ans-=sum[i]*(n-k);
            }
        printf("%.7lf\n",(double)ans/tot);
    }
    return 0;
}


评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值