「6月雅礼集训 2017 Day4」寻找天哥
【题目大意】
给出$n$个三维向量，设当前向量长度为$L$，每次沿着向量等概率走$[0,L]$个长度。一个球每秒半径增加1个长度，直到覆盖位置，每秒耗能为球体积，求总耗能的期望。
设最后半径为R，那么求得就是$ \int_0^R \frac{4}{3}\pi x^3\, dx.$的期望。
$1 \leq n \leq 3000$
【题解】
也就是求$E(\frac{\pi}{3}R^4)$，问题在于怎么求$E(R^4)$。
先提供一种错误做法及其实现：
我们设向量为$\{p_n\}$，设$x_i$是$(0,1)$等概率随机的。
那么相当于求$E( (\sum_{i=1}^n p_ix_i)^4 )$。
拆出一个数，相当于
$E((\sum_{i=1}^{n-1} p_ix_i + p_nx_n)^4)$
二项式展开，得
$E( (\sum_{i=1}^{n-1}p_ix_i)^4 + 4(\sum_{i=1}^{n-1}p_ix_i)^3(p_nx_n) + 6(\sum_{i=1}^{n-1}p_ix_i)^2(p_nx_n)^2+4(\sum_{i=1}^{n-1}p_ix_i)(p_nx_n)^3 + (p_nx_n)^4)$
所以我们只要维护$p_ix_i$的1~4次方的期望就行了吗？
讲道理是的啊，但是这种做法是错的，只能在所有向量同向的时候是对的。
为什么呢？因为考虑期望中有向量和向量数乘的一项，比如$4(\sum_{i=1}^{n-1}p_ix_i)(p_nx_n)^3$，这是不支持结合律的！！！
所以。。是错的qwq
======================分割线======================
我们接下来讲正确的做法：
我们考虑把向量分成三个坐标表示$(a_i,b_i,c_i)$。
求$E(R^4)$还可以看做$E( ((\sum_{i=1}^{n} a_ix_i)^2 + (\sum_{i=1}^{n} b_ix_i)^2 + (\sum_{i=1}^{n} c_ix_i)^2) ^2)$
这样好像正常多了，至少没有向量了。
设当前做到k，当前的$p = \sum_{i=1}^k a_ix_i$，$q = \sum_{i=1}^k b_ix_i$，$r = \sum_{i=1}^k c_ix_i$。
那么也就是求$E( (p^2+q^2+r^2)^2 ) = E(p^4+q^4+r^4+2p^2q^2+2p^2r^2+2q^2r^2)$
设f[x,i,j,k]表示加到第x个向量，$p^i * q^j * r^k$的期望。
那么根据期望的线性性，答案就是f[n,4,0,0]+f[n,0,4,0]+f[n,0,0,4]+2 * (f[n,0,2,2]+f[n,2,0,2]+f[n,2,2,0])
考虑转移，每次加入一个向量，我们试着把其中一项加入当前的$a_ix_i,b_ix_i,c_ix_i$（记为$p',q',r'$）。
$E((p+p')^2(q+q')^2) = E((p^2+2pp'+p'^2)(q^2+2qq'+q'^2)) = E(p^2q^2+ 2p^2qq' + p^2q'^2+2pq^2p' + 4pqp'q' + 2pp'q'^2 + q^2p'^2 + 2qp'^2q' + p'^2q'^2)$
可能已经发现了转移方式了
f[x,i,j,k] = f[x-1, i-A, j-B, k-C] * E(A, B, C) * C(i, A) * C(j, B) * C(k, C)
E(a,b,c)表示$p'^A * q'^B * r'^C$的期望，根据定义显然就是求$E(a_i^Ab_i^Bc_i^Cx_i^{A+B+C})$的期望。
这里$x_i$是$(0,1)$的变量，所以$E(x_i^p) = 1/(p+1)$。由于$x_i$和前面几个x都是互相独立的，所以这时候$E(AB) = E(A)E(B)$.
$a_i^Ab_i^Bc_i^C$都是常数，所以最后答案是$a_i^Ab_i^Bc_i^C / (A+B+C+1)$
然后就可以转移啦。。
复杂度O(n * 常数)
# include <math.h> # include <stdio.h> # include <string.h> # include <iostream> # include <algorithm> using namespace std; typedef long long ll; typedef unsigned long long ull; typedef long double ld; const int M = 3e3 + 10; const int mod = 1e9 + 7; const ld pi = acos(-1.0); int C[5][5], n; ld f[2][5][5][5], E[5][5][5]; ld a[M], b[M], c[M], sa[M], sb[M], sc[M]; // E(R^4) = E( ((∑a_i*x_i)^2 + (∑b_i*x_i)^2 + (∑c_i*x_i)^2)^2 ) // p = ∑a_i*x_i, q = ∑b_i*x_i, r = ∑c_i*x_i // E(R^4) = E( p^4 + q^4 + r^4 + 2p^2 q^2 + 2p^2 r^2 + 2q^2 r^2 ) // p = p + a_nx_n, q = q + b_nx_n, r = r + c_nx_n // E((a_nx_n) ^ A * (b_nx_n) ^ B * (c_nx_n) ^ C) /* e.g. (p+p')^2(q+q'^2) = (p^2+2pp'+p'^2)(q^2+2qq'+q'^2) = p^2q^2 + 2p^2q * q' + p^2q'^2 + ... */ int main() { // freopen("find.in", "r", stdin); // freopen("find.out", "w", stdout); C[0][0] = 1; for (int i=1; i<=4; ++i) { C[i][0] = 1; for (int j=1; j<=i; ++j) C[i][j] = C[i-1][j] + C[i-1][j-1]; } while(cin >> n) { if(!n) continue; memset(f, 0, sizeof f); int cur = 1, pre = 0; f[pre][0][0][0] = 1; double Alpha, Beta, L; for (int i=1; i<=n; ++i) { scanf("%lf%lf%lf", &Alpha, &Beta, &L); a[i] = L * cos(Beta) * cos(Alpha), b[i] = L * cos(Beta) * sin(Alpha), c[i] = L * sin(Beta); } for (int i=1; i<=n; ++i) { sa[0] = sb[0] = sc[0] = 1; for (int j=1; j<=4; ++j) { sa[j] = sa[j-1] * a[i]; sb[j] = sb[j-1] * b[i]; sc[j] = sc[j-1] * c[i]; } for (int i=0; i<=4; ++i) for (int j=0; i+j<=4; ++j) for (int k=0; i+j+k<=4; ++k) E[i][j][k] = sa[i] * sb[j] * sc[k] / (i+j+k+1); for (int i=0; i<=4; ++i) for (int j=0; i+j<=4; ++j) for (int k=0; i+j+k<=4; ++k) { f[cur][i][j][k] = 0; for (int x=0; x<=i; ++x) for (int y=0; y<=j; ++y) for (int z=0; z<=k; ++z) f[cur][i][j][k] += f[pre][i-x][j-y][k-z] * E[x][y][z] * C[i][x] * C[j][y] * C[k][z]; } swap(pre, cur); } ld ans = f[pre][4][0][0] + f[pre][0][4][0] + f[pre][0][0][4] + 2.0 * (f[pre][2][2][0] + f[pre][2][0][2] + f[pre][0][2][2]); ans = ans / 3.0 * pi; printf("%.9lf\n", (double)ans); } return 0; }
下面这份是只能过共线的代码：
# include <math.h> # include <stdio.h> # include <string.h> # include <iostream> # include <algorithm> using namespace std; typedef long long ll; typedef unsigned long long ull; typedef long double ld; const int M = 3e3 + 10, N = 1e5 + 10; const int mod = 1e9 + 7; const double pi = acos(-1.0); int n; struct vec { bool f; double a, b, c; vec() {} inline void set(bool _f, double _a, double _b = 0.0, double _c = 0.0) { f = _f; a = _a, b = _b, c = _c; } friend vec operator + (vec a, vec b) { vec ret; if(a.f && b.f) ret.set(1, a.a+b.a, a.b+b.b, a.c+b.c); else ret.set(0, a.a+b.a); return ret; } friend vec operator * (vec a, vec b) { vec ret; if(a.f && b.f) ret.set(0, a.a * b.a + a.b * b.b + a.c * b.c); if(a.f && !b.f) ret.set(1, a.a * b.a, a.b * b.a, a.c * b.a); if(!a.f && b.f) ret.set(1, a.a * b.a, a.a * b.b, a.a * b.c); if(!a.f && !b.f) ret.set(0, a.a * b.a); return ret; } friend vec operator * (vec a, double b) { vec ret; if(a.f) ret.set(1, a.a*b, a.b*b, a.c*b); else ret.set(0, a.a*b); return ret; } friend vec operator / (vec a, double b) { vec ret; if(a.f) ret.set(1, a.a/b, a.b/b, a.c/b); else ret.set(0, a.a/b); return ret; } inline void out() { if(f) printf("%lf, %lf, %lf\n", a, b, c); else printf("%lf\n", a); } }p[M][5], a[M], f[M][5]; int main() { freopen("find.in", "r", stdin); freopen("find.out", "w", stdout); while(cin >> n) { double Alpha, Beta, L; if(n == 0) break; for (int i=1; i<=n; ++i) { scanf("%lf%lf%lf", &Alpha, &Beta, &L); a[i].set(1, L * cos(Beta) * cos(Alpha), L * cos(Beta) * sin(Alpha), L * sin(Beta)); p[i][1] = a[i] / 2.0; p[i][2] = (a[i] * a[i]) / 3.0; p[i][3] = (a[i] * a[i] * a[i]) / 4.0; p[i][4] = (a[i] * a[i] * a[i] * a[i]) / 5.0; // cout << "i = " << i << endl; // p[i][1].out(); // p[i][2].out(); // p[i][3].out(); // p[i][4].out(); // cout << endl; } /* ans = 1/3 * pi * E(r^4) E((∑p[i]x[i])^4) = E((∑p[i]x[i] + p[n]x[n]) ^ 4) = E((∑p[i]x[i]) ^ 4) + 3E((∑p[i]x[i])^3)E(p[n]x[n]) + 6E((∑p[i]x[i])^2)E((p[n]x[n])^2) + 3E(∑p[i]x[i]) E((p[n]x[n])^3) + E((p[n]x[n])^4) E( (p[i]x[i])^3 ) = */ f[1][1] = p[1][1], f[1][2] = p[1][2], f[1][3] = p[1][3], f[1][4] = p[1][4]; for (int i=2; i<=n; ++i) { f[i][1] = f[i-1][1] + p[i][1]; f[i][2] = f[i-1][2] + (f[i-1][1] * p[i][1]) * 2 + p[i][2]; f[i][3] = f[i-1][3] + (f[i-1][2] * p[i][1]) * 3 + (f[i-1][1] * p[i][2]) * 3 + p[i][3]; f[i][4] = f[i-1][4] + (f[i-1][3] * p[i][1]) * 4 + (f[i-1][2] * p[i][2]) * 6 + (f[i-1][1] * p[i][3]) * 4 + p[i][4]; } double ans = f[n][4].a * pi / 3.0; printf("%.10lf\n", ans); } return 0; }
