【GDOI2018模拟7.7】寻找天哥

Description

n<=3000

Solution

答案相当于求

E(∫R04π3x3dx)

也就是求

E(π3R4)

现在就是要找到一种方法求E(R4)
发现这个东西等价于求

E(((∑i=1naixi)2+(∑i=1nbixi)2+(∑i=1ncixi)2)2)

其中xi是一个在[0,1]范围内均匀分布的实数
设∑ni=1aixi=A,∑ni=1bixi=B,∑ni=1cixi=C
那么答案就是E(A4+B4+C4+2(A2B2+A2C2+B2C2))
设Fx,i,j,k表示当前做到第x个向量，形如Ai∗Bj∗Ck的式子的期望
考虑新加入一个向量对这条式子的影响,设新加入的向量为(A’,B’,C’)
我们可以知道答案是由很多个形如这样的式子组成的

E((A+A′)2(B+B′)2)

那么根据二项式定理我们可以知道从(Ai,Bj,Ck)转移到(AiA′i′,BjB′j′,CkC′k′)是乘上一个Cii+i′Cjj+j′Ckk+k′的系数的
这样子我们就只需要对每个向量计算出E(aixbjxckxxi+j+k)，就可以很简便的转移了
显然∫10xpdx=1p+1
所以E(aixbjxckxxi+j+k)=aixbjxckxi+j+k+1
那么就直接Dp转移就好了，复杂度O(n*常数)

Code

#include <cmath>#include <cstdio>#include <cstring>#include <algorithm>#define fo(i,a,b) for(int i=a;i<=b;i++)using namespace std;typedef double db;const int N=3*1e3+5;int CC[5][5],n;db alpha,beta,l,a[N],b[N],c[N],f[N][5][5][5],e[5][5][5],sa[5],sb[5],sc[5];int main() {    freopen("tg.in","r",stdin);    fo(i,0,4) {        CC[i][0]=1;        fo(j,1,i) CC[i][j]=CC[i-1][j-1]+CC[i-1][j];    }    for(scanf("%d",&n);n;scanf("%d",&n)) {        fo(i,1,n) {            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);        }        memset(f,0,sizeof(f));        f[0][0][0][0]=1;        fo(x,1,n) {            sa[0]=sb[0]=sc[0]=1;            fo(i,1,4) {                sa[i]=sa[i-1]*a[x];                sb[i]=sb[i-1]*b[x];                sc[i]=sc[i-1]*c[x];            }            fo(i,0,4)                fo(j,0,4)                    fo(k,0,4)                         e[i][j][k]=sa[i]*sb[j]*sc[k]/(i+j+k+1);            fo(i,0,4)                fo(j,0,4)                    fo(k,0,4)                         fo(A,0,i)                            fo(B,0,j)                                fo(C,0,k)                                     f[x][i][j][k]+=f[x-1][i-A][j-B][k-C]*e[A][B][C]*CC[i][A]*CC[j][B]*CC[k][C];             }        db ans=f[n][4][0][0]+f[n][0][4][0]+f[n][0][0][4]+2*(f[n][2][2][0]+f[n][2][0][2]+f[n][0][2][2]);        ans=ans*M_PI/3.0;        printf("%.6lf\n",ans);    }}