Kattis - heapsoffun Heaps of Fun (概率密度函数+dp)

Posted 2022-08-14 asdfsag

题意：有一棵含有n个结点(n<=300)的根树，树上每个结点上的权值是从[0,ai](ai<=1e9)区间内随机的一个实数，问这棵树能形成一个最小堆的概率。

由于结点取值范围是1e9而且是实数，所以枚举权值dp自然是行不通的了，但可以从函数的角度上考虑。

首先需要了解两个概念：

CDF:分布函数，记为F(x)，表示函数F的取值小于等于x的概率。

PDF:概率密度函数，记为f(x)，是F(x)的导数，反之，F(x)是f(x)在区间(-∞,x]上的积分。由于本题所有的取值都是从0开始的，因此也可以表示在区间[0,x]上的积分。

那么答案就是根节点rt所对应的CDF的上界Frt(+∞)。

我们记fu(x)为结点u的取值的PDF，Fu(x)为其CDF，那么对于每个叶子结点u来说：

$f_u(x)=\left\\beginmatrix\beginaligned&\frac1a[u],0\leqslant x\leqslant a[u]\\&0,others\endaligned\endmatrix\right.,F_u(x)=\int_0^xf_u(t)dt=\left\\beginmatrix\beginaligned&\frac1a[u]x,0\leqslant x\leqslant a[u]\\&1,x>a[u]\endaligned\endmatrix\right.$

那么有了子结点的PDF和CDF，如何去求父结点的PDF和CDF呢？

首先我们要算出每个子结点的取值大于x的概率，因为如果结点u的某个取值x是合法的，那么它的所有子节点的取值都必须大于x。既然Fu(x)表示结点u取值小于等于x的概率，那么取值大于x的概率呢？是1-Fu(x)吗？No！这个公式只适用于叶子结点，因为这里的Fu(x)表示的是“该结点所在子树满足最小堆性质且该结点的值小于等于x的概率”，因此Fu的上界不一定为1，也就是说x的所有取值的概率之和不一定为1。那么应该如何计算呢？

设ub[u]表示结点u及子树下的所有结点的最小的a，那么显然结点u的取值不能超过ub[u]，即ub[u]是结点u取值的上界。

所以，正确的计算方法是结点u取值大于x的概率$G_u(x)=\left\\beginmatrix\beginaligned&F(ub[u])-F(x),0\leqslant x\leqslant ub[u]\\&0,x>ub[u]\endaligned\endmatrix\right.$

有了子结点的G，就可以求父结点的f了，通过概率计算公式可得：

$f_u(x)=\frac1a[u]\prod\limits_fa[v]=uF‘_v(x)$

然后对fu(x)积分可以得到Fu(x)，然后又可以求出Gu(x)，之后就可以继续往父结点递推了。最终答案就是根节点的G(x)的常数项。

代码如下：（为了方便，代码中每个结点最后保存的函数是G(x)）

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 typedef long long ll;
 4 const int N=300+10,mod=1e9+7;
 5 int hd[N],ne,rt,a[N],n,ub[N];
 6 struct E int v,nxt; e[N];
 7 void addedge(int u,int v) e[ne]= v,hd[u],hd[u]=ne++;
 8 int Pow(int x,int p) 
 9     int ret=1;
10     for(; p; p>>=1,x=(ll)x*x%mod)if(p&1)ret=(ll)ret*x%mod;
11     return ret;
12 
13 int inv(int x) return Pow(x,mod-2);
14 typedef vector<int> Poly;
15 Poly F[N];
16 Poly operator*(Poly& a,Poly& b) 
17     Poly c(a.size()+b.size()-1,0);
18     for(int i=0; i<a.size(); ++i)
19         for(int j=0; j<b.size(); ++j)
20             c[i+j]=(c[i+j]+(ll)a[i]*b[j])%mod;
21     return c;
22 
23 Poly itg(Poly& a) 
24     Poly c(a.size()+1,0);
25     for(int i=0; i<a.size(); ++i)c[i+1]=(ll)a[i]*inv(i+1)%mod;
26     return c;
27 
28 int eval(Poly& a,int x) 
29     int ret=0;
30     for(int i=a.size()-1; i>=0; --i)ret=((ll)ret*x+a[i])%mod;
31     return ret;
32 
33 void dfs(int u) 
34     ub[u]=a[u],F[u].push_back(inv(a[u]));
35     for(int i=hd[u]; ~i; i=e[i].nxt) 
36         int v=e[i].v;
37         dfs(v),ub[u]=min(ub[u],ub[v]),F[u]=F[u]*F[v];
38     
39     F[u]=itg(F[u]),F[u][0]=eval(F[u],ub[u]);
40     for(int i=1; i<F[u].size(); ++i)F[u][i]=-F[u][i];
41 
42 int main() 
43     memset(hd,-1,sizeof hd),ne=0;
44     scanf("%d",&n);
45     for(int i=1,f; i<=n; ++i) 
46         scanf("%d%d",&a[i],&f);
47         if(f)addedge(f,i);
48         else rt=i;
49     
50     dfs(rt);
51     printf("%d\n",(F[rt][0]+mod)%mod);
52     return 0;
53