【JZOJ100208】【20190705】传说之下

题目

三维空间上有一个点，进行了(n)次移动

第(i)次为在([0,L_i])内随机一个长度(l_i)，向(vec P_i)方向移动(l_i)

$vec P_i $ 表示为 ((alpha_i,eta_i)) ，意义为设 (vec P_i) 在 (xy) 上的投影为 (vec Q_i) ， (alpha_i) 为 (vec Q_i) 和 (xy) 的夹角，(eta_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) 　

  根据二项式定理

  [egin{align} &dp_{i,j,k,l} = sum_psum_qsum_r dp_{i-1,j-p,k-q,l-r} imes (^j_p)(^k_q)(^l_r) imes 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
",ans);
	}
	return 0;
}