【JZOJ100208】【20190705】傳說之下

題目

三維空間上有一個點，進行了\(n\)次移動

第\(i\)次爲在\([0,L_i]\)內隨機一個長度\(l_i\)，向\(\vec P_i\)方向移動\(l_i\)

$\vec P_i $ 表示爲 \((\alpha_i,\beta_i)\) ，意義爲設 \(\vec P_i\) 在 \(xy\) 上的投影爲 \(\vec Q_i\) ， \(\alpha_i\) 爲 \(\vec Q_i\) 和 \(xy\) 的夾角，\(\beta_i\) 爲 \(\vec P_i\) 和 $ \vec Q_i$ 的夾角

從原點開始有一個球，每秒半徑增長1個單位，在時刻\(i\)會消耗當前體積的能量

求消耗能量的指望值

$ n \le 3000 $

題解

• 考慮最後的半徑\(R\)，答案即 $ E (\int_0^R \frac{4}{3} \pi x^3 dx) = \frac{\pi}{3} E(R^4)$

• 設第\(i\)次移動的向量爲\((a_i,b_i,c_i)x_i\) ， \(x_i\) 爲一個在 \([0,1]\)隨機分佈的變量

• $E(R^4) =E ( ( (\sum a_ix_i)^2 + (\sum b_ix_i)^2 + (\sum c_ix_i)^2) ^2 ) $

• 設 $ A_i = \sum a_ix_i $ ,BC同理，設 $ dp_{i,j,k,l} = E(A_i^jB_i^kC_i^l) $

• 只須要求出\(dp\)便可求出 \(ans\) 　

  根據二項式定理
  \[ \begin{align} &dp_{i,j,k,l} = \sum_p\sum_q\sum_r dp_{i-1,j-p,k-q,l-r} \times (^j_p)(^k_q)(^l_r) \times a_i^pb_i^qc_i^rE(x_i^{p+q+r})\\ &因爲f(x)在[L,R]內的指望E(f(x)) = \frac{\int_L^R f(x) dx}{R-L} \\ &因此E(x_i^{p+q+r}) \ = \ \frac{1}{p+q+r+1} \\ \end{align} \]

• 時間複雜度：\(O(5^6n)\)

• 因爲有大量無用狀態，採用記憶化搜索

Code

#include<bits/stdc++.h>
#define ld long double
using namespace std;
const int N=3010;
int T,n,Tim,C[5][5],g[N][5][5][5];
ld a[N][5],b[N][5],c[N][5],f[N][5][5][5],ny[13];
ld cal(int i,int j,int k,int l){
    if(!i)return !j&&!k&&!l;
    if(g[i][j][k][l]==Tim)return f[i][j][k][l];
    g[i][j][k][l]=Tim;
    ld &res=f[i][j][k][l];
    res=0;
    for(int p=0;p<=j;++p)
        for(int q=0;q<=k;++q)
            for(int r=0;r<=l;++r)
                res+=cal(i-1,j-p,k-q,l-r)*C[j][p]*C[k][q]*C[l][r]*a[i][p]*b[i][q]*c[i][r]*ny[p+q+r+1];
    return res;
}
int main(){
    freopen("undertale.in","r",stdin);
    freopen("undertale.out","w",stdout);
    scanf("%d",&T);
    for(int i=0;i<5;++i)C[i][0]=1;
    for(int i=1;i<5;++i)
    for(int j=1;j<5;++j)
        C[i][j]=C[i-1][j]+C[i-1][j-1];
    for(int i=1;i<13;++i)ny[i]=1.0/i;
    while(T--){
        scanf("%d",&n);
        for(int i=1;i<=n;++i){
            ld x,y,l;
            scanf("%Lf%Lf%Lf",&x,&y,&l);
            a[i][0]=b[i][0]=c[i][0]=1;
            c[i][1]=sin(y)*l;l*=cos(y);
            b[i][1]=sin(x)*l;
            a[i][1]=cos(x)*l;
            for(int j=2;j<5;++j){
                a[i][j]=a[i][j-1]*a[i][1];
                b[i][j]=b[i][j-1]*b[i][1];
                c[i][j]=c[i][j-1]*c[i][1];
            }
        }
        ++Tim;
        ld ans = cal(n,4,0,0) + cal(n,0,4,0) + cal(n,0,0,4)
                +2*cal(n,2,2,0) + 2*cal(n,2,0,2) + 2*cal(n,0,2,2);
        ans*=acos(-1)/3;
        printf("%.10Lf\n",ans);
    }
    return 0;
}