# BZOJ-3771. Triple

Posted by miskcoo on April 23, 2015

 $x_1$ $x_2$ $x_3$ 1 2 3

 价值 1 2 3 4 5 6 方案数 1 1 2 1 1 1

$A(x) = x^1 + x^2 + x^3$

$B(x) = x^2 + x^4 + x^6 \\ C(x) = x^3 + x^6 + x^9$

$\frac{A^3(x) - 3A(x)\cdot B(x) + 2C(x)}{6}$

$\frac{A^2(x)-B(x)}{2}$

$A(x)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
/* BZOJ-3771: Triple
*   母函数+FFT */
#include <cstdio>
#include <cmath>
#include <complex>
#include <algorithm>
#include <iostream>

typedef std::complex<double> complex_t;

void shuffle(int n, complex_t *x)
{
for(int i = 0, j = 0; i != n; ++i)
{
if(i > j) std::swap(x[i], x[j]);
for(int l = n >> 1; (j ^= l) < l; l >>= 1);
}
}

void transform(int n, complex_t *x, complex_t *w)
{
shuffle(n, x);
for(int i = 2; i <= n; i <<= 1, ++w)
{
int m = i >> 1;
complex_t *tmp = new complex_t[m];
tmp[0] = 1;
for(int p = 1; p != m; ++p)
tmp[p] = tmp[p - 1] * *w;

for(int p = 0; p < n; p += i)
{
for(int t = 0; t != m; ++t)
{
complex_t z = x[p + m + t] * tmp[t];
x[p + m + t] = x[p + t] - z;
x[p + t] += z;
}
}

delete[] tmp;
}
}

const int MaxL = 17, MaxN = (1 << MaxL) + 10;
complex_t eps[MaxL + 1], eps_inv[MaxL + 1];
complex_t A[MaxN], B[MaxN], C[MaxN];
void init_eps()
{
double pi = std::acos(-1.0);
double angle = 2.0 * pi / (1 << MaxL);
eps[MaxL] = complex_t(std::cos(angle), std::sin(angle));
for(int i = MaxL - 1; i >= 0; --i)
eps[i] = eps[i + 1] * eps[i + 1];
for(int i = 0; i <= MaxL; ++i)
eps_inv[i] = 1.0 / eps[i];
}

int main()
{
init_eps();
int n, max_v = 0;
std::scanf("%d", &n);
for(int i = 0; i != n; ++i)
{
int v;
std::scanf("%d", &v);
A[v] += 1.0;
B[v * 2] += 1.0;
C[v * 3] += 1.0;
if(v * 3 > max_v) max_v = v * 3;
}

int m = 1;
while(m < max_v) m <<= 1;
transform(m, A, eps + 1);
transform(m, B, eps + 1);
for(int i = 0; i != m; ++i)
A[i] = A[i] * A[i] * A[i] / 6.0 + (A[i] * (A[i] - B[i]) - B[i]) / 2.0 + A[i];
transform(m, A, eps_inv + 1);
for(int i = 0; i != m; ++i)
A[i] = A[i] / double(m) + C[i] / 3.0;
for(int i = 0; i != m; ++i)
{
int x = int(A[i].real() + 0.5);
if(x) std::printf("%d %d\n", i, x);
}
return 0;
}