2018-01-20 HDU 4609 FFT 快速傅里叶变换

终于来了! 这道题目.

题目描述:
King OMeGa catched three men who had been streaking in the street. Looking as idiots though, the three men insisted that it was a kind of performance art, and begged the king to free them. Out of hatred to the real idiots, the king wanted to check if they were lying. The three men were sent to the king's forest, and each of them was asked to pick a branch one after another. If the three branches they bring back can form a triangle, their math ability would save them. Otherwise, they would be sent into jail.
However, the three men were exactly idiots, and what they would do is only to pick the branches randomly. Certainly, they couldn't pick the same branch - but the one with the same length as another is available. Given the lengths of all branches in the forest, determine the probability that they would be saved.

输入:
An integer T(T≤100) will exist in the first line of input, indicating the number of test cases.
Each test case begins with the number of branches N(3≤N≤10⁵).
The following line contains N integers a_i (1≤a_i≤10⁵), which denotes the length of each branch, respectively.

输出:
Output the probability that their branches can form a triangle, in accuracy of 7 decimal places.

样例输入:
2
4
1 3 3 4
4
2 3 3 4

样例输出:
0.5000000
1.0000000

题目的意思很简单, 给你一堆长度各异的棍子, 棍子的数量和每根棍子有多长也告诉你, 问你随便拿三根可以组成三角形的概率是多少, (所以3个滞胀只是用来吐槽的, 并没有什么实际的题意)

如果暴力当然是超时啦...

再怎么优化都没用的,肯定超时....

除非你有什么黑科技.....

呃, 一个不那么黑科技的东西就是使用FFT来优化.

如何优化, 当然是用多项式了, 想一下生成函数, 可以获得一个任取两根棍子获得长度的生成数组

首先先是FFT来进行多项式相乘, 然后得到一个生成函数的系数数组.

但是, 里面会有重复的, 这里是第一个去重, 比较容易想到

先想一下这个多项式相乘什么意思, 就是两堆一模一样的棍子, 第一堆随便拿一根, 第二堆随便拿一根

那么会多出来这么两种情况 1.两堆里面拿出来的两根棍子是一样的 2.第一堆拿了第a根棍子 第二堆拿了第b根棍子 这种情况和 第一堆拿了第b根棍子 第二堆拿了第a根棍子 是一样的.

针对情况一, 先遍历这一堆里面的每一根棍子, 把长度*2的次数(生成函数的系数)减一, 解决了情况1的重复

然后, 把所有的项的系数都除以2, 就解决了情况2的重复.

然后我们就得到了一个数组, 就是在一堆里面任取两根长度为何的生成数组(times[5]=4的意思是 任取两根棍子长度恰为5的方案数是4)

不过说实话这只是万里长征第一步, 后面的去重才麻烦着呢......

由于前面的操作, 我们要对原始的存储棍子的数组进行排序. 这里对每根棍子进行遍历, 求针对每根棍子有几种(不重复的)方案.

遍历时, 针对特定的棍子(比如到第i根棍子), 先不妨假设这跟棍子是最长的.

当然了,不管它是不是最长的, 它的长度一定是短于(小于!)其他两根棍子长度之和的, 所以把生成函数系数数组里面下标大于这根特定的棍子长度的数据加起来, 这里为了减少时间, 我们先处理出前缀和数组sum[], 这样结果就是sum[len]-sum[data[i]], 其中len是最大脚标, data[i]是特定棍子的长度.

如何保证它是不重复的呢? 关键就是前面提到的"不妨假设这跟棍子是最长的"

那么保证去重的正确性就是保证这根棍子是最长的. 对于刚刚算出来的sum[len]-sum[data[i]], 里面有一些错误的方案, 这是要去掉的.(下面的式子都是保证已经对data数组排过序之后的从0~n-1之间的遍历!)

1.这跟棍子并不是最大的之---它是最小的, 这种错误的方案数是 (n-i-1)*(n-i-2)/2

2.这根棍子并不是最大的之---有一根比它长, 另一根比它短, 这种错误的方案数是 (n-i-1)*i

3.咦, 为什么拿出来的两根棍子有一根也是第i根?---组合里面会有包含自己的, 这种错误的方案数是 n-1

至于方案数为什么是这么些, 这就是你组合数学的功底了, 跟我没关系.

然后减掉就行.

还有就是, 关于强制类型转换, 我真的心累啊.....一个都不能不一样的啊......

代码如下:

#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <algorithm>
using namespace std;
const double PI=acos(-1);

struct complex
{
    double a,b;//a表示实部,b表示虚部
    complex(double aa=0.0,double bb=0.0)
    {
        a=aa;
        b=bb;
    }
    complex operator +(const complex &e)
    {
        return complex(a+e.a,b+e.b);
    }
    complex operator -(const complex &e)
    {
        return complex(a-e.a,b-e.b);
    }
    complex operator *(const complex &e)
    {
        return complex(a*e.a-b*e.b,a*e.b+b*e.a);
    }
};

void change(complex *y,int len)
{
    int i,j,k;
    for(i=1,j=len/2;i<len-1; i++)
    {
        if(i<j)
            swap(y[i],y[j]);
        k=len/2;
        while(j>=k)
        {
            j-=k;
            k>>=1;
        }
        if(j<k)
        {
            j+=k;
        }
    }
}
void fft(complex *y,int len,int on)
{
    change(y, len);
    int i,j,k;
    for(i=2;i<=len;i<<=1)
    {
        complex wn(cos(-on*2*PI/i),sin(-on*2*PI/i));
        for(j=0;j<len;j+=i)
        {
            complex w(1,0);
            for(k=j;k<j+i/2;k++)
            {
                complex u=y[k],t=w*y[k+i/2];
                y[k]=u+t;
                y[k+i/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1)
        for(i=0;i<len;i++)
            y[i].a/=len;
}

const int maxx=500000;
complex a[maxx],b[maxx];
long long int sum[maxx];
int data[maxx];
long long int times[maxx];
int len,len_a;
//怕函数体里面塞不下所以.....

int main()
{
    int T,n,i;
    double res;
    long long int yes,all;
    scanf("%d",&T);
    while(T--)
    {
        yes=all=res=0;
        len=1;
        len_a=0;
        memset(times,0,sizeof(times));
        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));
        scanf("%d",&n);
        for(i=0;i<n;i++)
        {
            scanf("%d",&data[i]);
            times[data[i]]++;
            len_a=max(len_a,data[i]);
        }
        len_a++;
        sort(data,data+n);
        while(len<(2*len_a))
            len=len<<1;
        for(i=0;i<len;i++)
            a[i].a=times[i];
        fft(a,len,1);
        for(i=0;i<len;i++)
            b[i]=a[i]*a[i];
        fft(b,len,-1);
        for(i=0;i<len;i++)
            times[i]=(long long int)(0.5+b[i].a);
        len=2*data[n-1];//更新len!
        for(i=0;i<n;i++)
            times[data[i]*2]--;
        for(i=0;i<=len;i++)//这里中间为什么是等于?
            times[i]=times[i]/2;
        sum[0]=times[0];
        for(i=1;i<=len;i++)//这里中间为什么也是等于?
            sum[i]=sum[i-1]+times[i];//求前缀和
        for(i=0;i<n;i++)
        {
            yes+=sum[len]-sum[data[i]];
            yes-=(long long)(n-1-i)*i;
            yes-=(n-1);
            yes-=(long long)(n-1-i)*(n-i-2)/2;//关键是这里的去重和去掉的不符合题目的部分!
        }
        all=(long long)n*(n-1)*(n-2)/6;
        res=(double)yes/all;
        printf("%.7lf\n",res);
    }
    return 0;
}