题目

三维空间上有一个点，进行了\(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;
}