【GDOI2018模拟7.7】寻找天哥

最新推荐文章于 2021-07-10 16:01:49 发布

alan_cty

最新推荐文章于 2021-07-10 16:01:49 发布

阅读量554

点赞数 1

分类专栏：计算几何概率与期望文章标签： GDOI2018模拟寻找天哥概率与期望积分计算几何

本文链接：https://blog.csdn.net/alan_cty/article/details/74779546

版权

概率与期望同时被 2 个专栏收录

18 篇文章 0 订阅

订阅专栏

计算几何

9 篇文章 0 订阅

订阅专栏

Description

这里写图片描述
n<=3000

Solution

答案相当于求

E (\int R 0 4 π 3 x 3 d x)

$E(\int_{0}^{R}{4\pi \over 3}x^3dx)$

也就是求

E (π 3 R 4)

$E({\pi \over 3}R^4)$
现在就是要找到一种方法求

E(R4) $E(R^4)$
发现这个东西等价于求

E (((\sum i = 1 n a i x i) 2 + (\sum i = 1 n b i x i) 2 + (\sum i = 1 n c i x i) 2) 2)

$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)$
其中xi是一个在[0,1]范围内均匀分布的实数
设

∑ni=1aixi=A $\sum_{i=1}^{n}a_ix_i=A$ ,

∑ni=1bixi=B $\sum_{i=1}^{n}b_ix_i=B$ ,

∑ni=1cixi=C $\sum_{i=1}^{n}c_ix_i=C$
那么答案就是

E(A4+B4+C4+2(A2B2+A2C2+B2C2)) $E(A^4+B^4+C^4+2(A^2B^2+A^2C^2+B^2C^2))$
设

Fx,i,j,k $F_{x,i,j,k}$ 表示当前做到第x个向量，形如

Ai∗Bj∗Ck $A^i*B^j*C^k$ 的式子的期望
考虑新加入一个向量对这条式子的影响,设新加入的向量为(A’,B’,C’)
我们可以知道答案是由很多个形如这样的式子组成的

E ((A + A') 2 (B + B') 2)

$E((A+A')^2(B+B')^2)$
那么根据二项式定理我们可以知道从

(Ai,Bj,Ck) $(A^i,B^j,C^k)$ 转移到

(AiA′i′,BjB′j′,CkC′k′) $(A^iA'^{i'},B^jB'^{j'},C^kC'^{k'})$ 是乘上一个

Cii+i′Cjj+j′Ckk+k′ $C_{i+i'}^{i}C_{j+j'}^{j}C_{k+k'}^{k}$ 的系数的
这样子我们就只需要对每个向量计算出

E(aixbjxckxxi+j+k) $E(a_{x}^ib_{x}^jc_{x}^kx^{i+j+k})$ ，就可以很简便的转移了
显然

∫10xpdx=1p+1 $\int_{0}^{1}x^pdx={1\over {p+1}}$
所以

E(aixbjxckxxi+j+k)=aixbjxckxi+j+k+1 $E(a_{x}^ib_{x}^jc_{x}^kx^{i+j+k})={a_{x}^ib_{x}^jc_{x}^k\over {i+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);
    }
}