「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
", (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
", a, b, c);
        else printf("%lf
", 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
", ans);
    }
    
    return 0;
}