【FFT】HDU 4609 3-idiots

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

题意：n条边长，求任选3条能组成三角形的概率。

思路：求出res后，减去一大一下，减去都大于它的边，减去该边和其他边组成的方案。

代码：

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <algorithm>
  5
  6 const int N =400007;
  7
  8 typedef long long ll;
  9
 10 struct comp {
 11     double x, y;
 12     comp(double _x=0, double _y=0) : x(_x), y(_y) {}
 13 };
 14
 15 namespace FFT {
 16     const double pi2 = acos(-1.) * 2;
 17     comp a[N], b[N], tmp[N];
 18     int n, bn;
 19     ll res[N];
 20     inline comp W(int n, bool inv) {
 21         double ang = inv ? -pi2 / n : pi2 / n;
 22         return comp(cos(ang), sin(ang));
 23     }
 24     int bitrev(int x) {
 25         int ans = 0;
 26         for (int i=1; i<=bn; ++i)
 27             ans <<= 1, ans |= x & 1, x >>= 1;
 28         return ans;
 29     }
 30     void dft(comp *a,bool inv) {
 31         int step, to; comp w, wi, A, B;
 32         for (int i=0; i<n; ++i) {
 33             to = bitrev(i);
 34             if (to > i) std::swap(a[to], a[i]);
 35         }
 36         for (int i=1; i<=bn; ++i) {
 37             wi = W(1<<i, inv); w = comp(1, 0);
 38             step = 1 << (i-1);
 39             for (int k=0; k<step; ++k) {
 40                 for (int j=0; j<n; j+=1<<i) {
 41                     int t = j | k, d = j|k|step;
 42                     A = a[t];
 43                     B.x  = w.x * a[d].x - w.y * a[d].y;
 44                     B.y  = w.x * a[d].y + w.y * a[d].x;
 45                     a[t].x = A.x + B.x, a[t].y = A.y + B.y;
 46                     a[d].x = A.x - B.x, a[d].y = A.y - B.y;
 47                 }
 48                 comp tmp;
 49                 tmp.x = w.x * wi.x - w.y * wi.y;
 50                 tmp.y = w.x * wi.y + w.y * wi.x;
 51                 w = tmp;
 52             }
 53         }
 54     }
 55     int mul(int n1, int *x1, int n2, int *x2) {
 56         n = std::max(n1, n2);
 57         for (bn = 0; (1<<bn) < n; ++bn); ++bn;
 58         n = 1 << bn;
 59         for (int i=0; i<n; ++i) a[i] = b[i] = comp(0, 0);
 60         for (int i=0; i<n1; ++i) a[i] = comp(x1[i], 0);
 61         for (int i=0; i<n2; ++i) b[i] = comp(x2[i], 0);
 62         dft(a, false); dft(b, false);
 63         for (int i=0; i<n; ++i) {
 64             tmp[i].x = a[i].x * b[i].x - a[i].y * b[i].y;
 65             tmp[i].y = a[i].x * b[i].y + a[i].y * b[i].x;
 66         }
 67         dft(tmp, true);
 68         for (int i=0; i<n; ++i) res[i] = (ll)(tmp[i].x/n + 0.1);
 69         for (--n; n && !res[n]; --n);
 70         return n+1;
 71     }
 72 };
 73
 74 int n;
 75 int a[N], num[N];
 76 ll sum[N];
 77
 78 int main() {
 79     int T;
 80     scanf("%d", &T);
 81     while (T-- > 0) {
 82         scanf("%d", &n);
 83         int up = 0;
 84         for (int i = 0; i < n; ++i) scanf("%d", &a[i]), up = std::max(up, a[i]);
 85         memset(num, 0, sizeof num);
 86         for (int i = 0; i < n; ++i) ++num[a[i]];
 87         int len = FFT::mul(up + 1, num, up + 1, num);
 88         for (int i = 0; i < n; ++i) {
 89             --FFT::res[a[i] + a[i]];
 90         }
 91         for (int i = 0; i < len; ++i) {
 92             FFT::res[i] >>= 1;
 93         }
 94         sum[0] = FFT::res[0];
 95         for (int i = 1; i < len; ++i) {
 96             sum[i] = sum[i - 1] + FFT::res[i];
 97         }
 98         ll ans = 0;
 99         for (int i = 0; i < n; ++i) {
100             ans += sum[len - 1] - sum[a[i]];
101             ans -= n - 1, ans -= (ll) i * (n - 1 - i);
102             ans -= (ll) (n - 1 - i) * (n - 2 - i) / 2;
103         }
104         ll tot = (ll) n * (n - 1) * (n - 2) / 6;
105         printf("%.7f\n", 1. * ans / tot);
106     }
107     return 0;
108 }