3-idiots hdu4609 母函数+FFT 组合数学题

时间：2023-12-16 16:20:14

http://acm.hdu.edu.cn/showproblem.php?pid=4609

题意：1e5个数，求取三个数能形成三角形的概率。

题解（这怎么会是fft入门题QAQ）：

　　概率的算法就是三角形取法/总取法。总取法就是C（n，3）.

　　三角形取法如何计算？

part1：构造母函数F（日常套路），每一项的次数为长度，系数为该长度的木棍数量，用FFT算F^2 ，

　　　　得到的多项式就包含了任意取两跟棍子得到的所有长度的方案数：其中次数为两根棍长之和，系数为该长度的方案数，

part2：去重，考虑part1中得到的系数，并非方案数（举个例子，（a+b）^2==a^2+b^2+2ab）

　　　　首先，对于每个平方项，因为棍子不能重复使用，所以对于k^m要减去平方前k^(m/2）的系数，

　　　　其次，对于剩下的每一项，由于多项式乘法每个乘法都做了两遍，所以得到的系数为方案数的两倍。要/2.

　　　　我们将处理过的系数存入数组A。

part3：算方案数：为了避免重复，对于每一根的木棍，计算以它为最长边的三角形方案数：为该长度到最大长度的A的系数和，（前缀和O（1）算出）

　　　　

　　　　减去其中包含它的方案数，为1*（n-1）.

　　　　减去其中它不是最长的方案数，为（n-1-1）*比它长的木棍数-两根都比它长的方案数/2（乘法原理多算了一遍，容斥减去），也可以理解为 一根比它长一根比他短+两根都比他长/2.（其中比它长的木棍数并不需要树状数组，只要排个序，他就是n-i）

end

总之，fft只是一个开始，后面巨麻烦orz

坑：套miskcoo dalao的fft，结果发现别人的是魔改版的，需要深刻地理解fft才会用orz

ac代码：

#include<iostream>

#include<sstream>

#include<fstream>

#include<vector>

#include<list>

#include<deque>

#include<queue>

#include<stack>

#include<map>

#include<set>

#include<bitset>

#include<algorithm>

#include<cstdio>

#include<cstdlib>

#include<cstring>

#include<cctype>

#include<cmath>

#include<ctime>

#include<iomanip>

#include<complex>

using namespace std;

//const double eps(1e-8);

typedef long long lint;

const double PI = acos(-1.0);

const int MaxL = , MaxN =  << MaxL;

typedef complex<double> complex_t;

complex_t f[MaxN], g[MaxN];

complex_t eps[MaxN], inv_eps[MaxN];

void init_eps(int p)

{

    double pi = acos(-);

    //double angle = 2.0 * pi / p;

    for (int i = ; i != p; ++i)

        eps[i] = complex_t(cos(2.0 * pi * i / p), sin(2.0 * pi * i / p));

    for (int i = ; i != p; ++i)

        inv_eps[i] = conj(eps[i]);

}

void transform(int n, complex_t *x, complex_t *w)

{

    for (int i = , j = ; i != n; ++i)

    {

        if (i > j) std::swap(x[i], x[j]);

        for (int l = n >> ; (j ^= l) < l; l >>= );

    }

    for (int i = ; i <= n; i <<= )

    {

        int m = i >> ;

        for (int j = ; j < n; j += i)

        {

            for (int k = ; k != m; ++k)

            {

                complex_t z = x[j + m + k] * w[n / i * k];

                x[j + m + k] = x[j + k] - z;

                x[j + k] += z;

            }

        }

    }

}

int branch[];

int num[];

complex_t a[];//(1 << 17 = 131072, 1 << 18 = 262144)

lint A[];

lint sumA[];//表示A[i]的前缀和

int main()

{

    int T, n;

    scanf("%d", &T);

    while (T--)

    {

        scanf("%d", &n);

        int maxBranch = ;

        for (int i = ; i < n; i++)

        {

            scanf("%d", branch + i);

            maxBranch = max(maxBranch, branch[i]);

        }

        memset(num, , sizeof(num));

        for (int i = ; i < n; i++)

            num[branch[i]]++;

        for (int i = ; i <= maxBranch; i++)

            a[i] = num[i];

        int len = ;

        while (len <= maxBranch) len <<= ;

        len <<= ;

        for (int i = maxBranch + ; i < len; i++)

            a[i] = 0.0;

        init_eps(len);

        transform( len,a ,eps);

        for (int i = ; i < len; i++)

            a[i] = a[i] * a[i];

        transform( len,a, inv_eps);

        for (int i = ; i <=  * maxBranch; i++)

            A[i] = (lint)(a[i].real() + 0.5)/len;

        for (int i = ; i <=  * maxBranch; i += )

            A[i] -= num[i >> ];

        for (int i = ; i <=  * maxBranch; i++)

            A[i] /= ;

        //到现在为止A[i]表示的是取两根不同的branch的长度和为i的组合种数

        sumA[] = ;

        for (int i = ; i <=  * maxBranch; i++)

            sumA[i] = sumA[i - ] + A[i];

        lint ans = ;

        sort(branch, branch + n);

        for (int i = ; i <= n; i++)//以第i根作为边最长的

        {

            lint tmp = sumA[ * maxBranch] - sumA[branch[i]];//另外两条边长度和要大于branch[i]

            tmp -= (lint)(n - i)*(n - );//比它长

            tmp += (lint)(n - i)*(n - i - ) / ;//两条都比他长

            tmp -= n - ;//另外两条的组合中包括它自己的组合

            ans += tmp;

        }

        double p = ans * . / n / (n - ) / (n - );

        printf("%.7f\n", p);

    }

    cin >> n;

    return ;

}

/*

*/

/*

2

4

1 3 3 4

4

2 3 3 4

*/

相关文章