HDU 4609 3-idiots（FFT）

传送门

3-idiots

Time Limit: 10000/5000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 5910    Accepted Submission(s): 2056

Problem Description

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.

 

Input

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

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

 

Sample Input

  
  

   
   2
4
1 3 3 4
4
2 3 3 4
  
  

 

Sample Output

  
  

   
   0.5000000
1.0000000
  
  

 

Source

2013 Multi-University Training Contest 1

 

Recommend

liuyiding   |   We have carefully selected several similar problems for you:   6170  6169  6168  6167  6166 

其实题目是给了n条线段。问随机取三个，可以组成三角形的概率。


 


其实就是要求n条线段，选3条组成三角形的选法有多少种。


 


首先题目给了a数组，


如样例一：


4


1 3 3 4


把这个数组转化成num数组，num[i]表示长度为i的有num[i]条。


样例一就是


num = {0   1   0    2    1}


代表长度0的有0根，长度为1的有1根，长度为2的有0根，长度为3的有两根，长度为4的有1根。


使用FFT解决的问题就是num数组和num数组卷积。


num数组和num数组卷积的解决，其实就是从{1 3 3 4}取一个数，从{1 3 3 4}再取一个数，他们的和每个值各有多少个


例如{0 1 0 2 1}*{0 1 0 2 1} 卷积的结果应该是{0 0  1  0  4  2  4  4  1 }


长度为n的数组和长度为m的数组卷积，结果是长度为n+m-1的数组。


 


{0 1 0 2 1}*{0 1 0 2 1} 卷积的结果应该是{0 0  1  0  4  2  4  4  1 }。


这个结果的意义如下：


从{1 3 3 4}取一个数，从{1 3 3 4}再取一个数


取两个数和为 2 的取法是一种：1+1


           和为 4 的取法有四种：1+3， 1+3  ，3+1 ，3+1


           和为 5 的取法有两种：1+4 ，4+1；


           和为 6的取法有四种：3+3,3+3,3+3,3+3,3+3


           和为 7 的取法有四种： 3+4,3+4,4+3,4+3


           和为 8 的取法有 一种：4+4


 


利用FFT可以快速求取循环卷积，具体求解过程不解释了，就是DFT和FFT的基本理论了。


总之FFT就是快速求到了num和num卷积的结果。只要长度满足>=n+m+1.那么就可以用循环卷积得到线性卷积了。

假设a[i]=L3
你想要的是L1<L2<L3的组合。
因为之前已经去过重了，所以不用考虑L1=L2的情况。
只需考虑，这下面三种情况

 1. L1<L3<L2
因为L3<L2，所以一定满足L1+L2>L3
这种情况离散化后大于a[i]数有n-i-1个，小于a[i]的数有i个(注意i从0开始)
所以cnt-=(n-1-i)*i;
2,  L1<L2=L3 或 L1=L3<L2
这种，情况下一条边只能取L3,而另一条边任意取都一定满足L1+L2>L
所以cnt-=n-1;

3  L3< L1<L2时
大于a[i]的数有n-i-1种
L1,L2从其中任选两种
这种情况下，因为L1一定要大于L2所以
总共有 （n-i-1）*（n-i-2）/2种
所以cnt-=（n-i-1）*（n-i-2）/2;

于是就有了

LL ans=0;
FOr(0,n,i)
{
    ans+=sum[len]-sum[a[i]];
    ans-=(LL)(n-i-1)*i;
    ans-=(n-1);
    ans-=(LL)(n-i-2)*(n-i-1)/2;
}

//china no.1
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <vector>
#include <iostream>
#include <string>
#include <map>
#include <stack>
#include <cstring>
#include <queue>
#include <list>
#include <stdio.h>
#include <set>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <cctype>
#include <sstream>
#include <functional>
#include <stdlib.h>
#include <time.h>
#include <bitset>
using namespace std;

#define pi acos(-1)
#define PI acos(-1)
#define endl '\n'
#define srand() srand(time(0));
#define me(x,y) memset(x,y,sizeof(x));
#define foreach(it,a) for(__typeof((a).begin()) it=(a).begin();it!=(a).end();it++)
#define close() ios::sync_with_stdio(0); cin.tie(0);
#define FOR(x,n,i) for(int i=x;i<=n;i++)
#define FOr(x,n,i) for(int i=x;i<n;i++)
#define W while
#define sgn(x) ((x) < 0 ? -1 : (x) > 0)
#define bug printf("***********\n");
#define db double
#define ll long long
typedef long long LL;
const int INF=0x3f3f3f3f;
const LL LINF=0x3f3f3f3f3f3f3f3fLL;
const int dx[]={-1,0,1,0,1,-1,-1,1};
const int dy[]={0,1,0,-1,-1,1,-1,1};
const int maxn=1e3+10;
const int maxx=4e5+100;
const double EPS=1e-8;
const double eps=1e-8;
const int mod=1e9+7;
template<class T>inline T min(T a,T b,T c) { return min(min(a,b),c);}
template<class T>inline T max(T a,T b,T c) { return max(max(a,b),c);}
template<class T>inline T min(T a,T b,T c,T d) { return min(min(a,b),min(c,d));}
template<class T>inline T max(T a,T b,T c,T d) { return max(max(a,b),max(c,d));}
template <class T>
inline bool scan_d(T &ret){char c;int sgn;if (c = getchar(), c == EOF){return 0;}
while (c != '-' && (c < '0' || c > '9')){c = getchar();}sgn = (c == '-') ? -1 : 1;ret = (c == '-') ? 0 : (c - '0');
while (c = getchar(), c >= '0' && c <= '9'){ret = ret * 10 + (c - '0');}ret *= sgn;return 1;}

inline bool scan_lf(double &num){char in;double Dec=0.1;bool IsN=false,IsD=false;in=getchar();if(in==EOF) return false;
while(in!='-'&&in!='.'&&(in<'0'||in>'9'))in=getchar();if(in=='-'){IsN=true;num=0;}else if(in=='.'){IsD=true;num=0;}
else num=in-'0';if(!IsD){while(in=getchar(),in>='0'&&in<='9'){num*=10;num+=in-'0';}}
if(in!='.'){if(IsN) num=-num;return true;}else{while(in=getchar(),in>='0'&&in<='9'){num+=Dec*(in-'0');Dec*=0.1;}}
if(IsN) num=-num;return true;}

void Out(LL a){if(a < 0) { putchar('-'); a = -a; }if(a >= 10) Out(a / 10);putchar(a % 10 + '0');}
void print(LL a){ Out(a),puts("");}
//freopen( "in.txt" , "r" , stdin );
//freopen( "data.txt" , "w" , stdout );
//cerr << "run time is " << clock() << endl;

//复数结构体
struct complex
{
    double r,i;
    complex(double _r = 0.0,double _i = 0.0)
    {
        r = _r; i = _i;
    }
    complex operator +(const complex &b)
    {
        return complex(r+b.r,i+b.i);
    }
    complex operator -(const complex &b)
    {
        return complex(r-b.r,i-b.i);
    }
    complex operator *(const complex &b)
    {
        return complex(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};
/*
 * 进行FFT和IFFT前的反转变换。
 * 位置i和 （i二进制反转后位置）互换
 * len必须去2的幂
 */
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]);
        //交换互为小标反转的元素，i<j保证交换一次
        //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的
        k = len/2;
        while( j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k) j += k;
    }
}
/*
 * 做FFT
 * len必须为2^k形式，
 * on==1时是DFT，on==-1时是IDFT
 */
void fft(complex y[],int len,int on)
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0;j < len;j+=h)
        {
            complex w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                complex u = y[k];
                complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0;i < len;i++)
            y[i].r /= len;
}
int a[maxx/4];
complex x1[maxx];
LL sum[maxx],num[maxx];
int main()
{

    int t;
    scan_d(t);
    W(t--)
    {
        int n;
        scan_d(n);
        me(num,0);
        FOr(0,n,i)
        {
            scan_d(a[i]);
            num[a[i]]++;
        }
        sort(a,a+n);
        int len1=a[n-1]+1;
        int len=1;
        W(len<len1*2)
            len<<=1;
        FOr(0,len1,i)
        {
            x1[i]=complex(num[i],0);
        }
        FOr(len1,len,i)
        {
            x1[i]=complex(0,0);
        }
        fft(x1,len,1);
        FOr(0,len,i)
        {
            x1[i]=x1[i]*x1[i];
        }
        fft(x1,len,-1);
        FOr(0,len,i)
        {
            num[i]=(LL)(x1[i].r+0.5);
        }
        //len=2*a[n-1];
        FOr(0,n,i)
        {
            num[a[i]+a[i]]--;
        }
        FOR(0,len,i)//
        {
            num[i]/=2;
        }
        sum[0]=0;
        FOR(1,len,i)
        {
            sum[i]=sum[i-1]+num[i];
        }
        LL ans=0;
        FOr(0,n,i)
        {
            ans+=sum[len]-sum[a[i]];
            ans-=(LL)(n-i-1)*i;
            ans-=(n-1);
            ans-=(LL)(n-i-2)*(n-i-1)/2;
        }
        LL tot=(LL)n*(n-1)*(n-2)/6;
        printf("%.7f\n",(db)ans/tot);
    }

}