【集训队互测2015】胡策的小树

概率 dp,树上高斯消元

题目大意:

给定一棵有根树,以随机父亲的方式随机生成,每个点有一个点权 $a_i$,保证 $a$ 是一个 $0\sim n-1$ 的排列。

有 $n$ 只猴子,初始每个结点各有一只。每秒钟,每个猴子都会进行移动。结点 $i$ 上的猴子有 $p_i=\frac {a_i}n$ 的概率会跳到它的父亲处(规定 $p_1=0$),如果没有成功跳到它父亲处,则会等概率随机跳到 $i$ 的子树中的一个结点(包括本身)。

令 $g_i$ 表示第 $i$ 秒钟时,成功跳到父亲结点的猴子占的比例,现在假设经过了无穷秒,求所有 $g$ 的平均值的期望。

现在要求选定一个 $x$($0\le x\lt n$),并令所有的 $a_i$ 变为 $(a_i+x)\bmod n$,问能得到的最大的期望是多少。

题解:

先考虑 $x=0$ 的情况。

由于轮数有无穷多,每只猴子都会以 1 的概率在有限的步数内跳到根,因此每只猴子对期望的贡献是一样的。

令 $x_i$ 表示某一时刻猴子在 $i$ 处的概率,则每只猴子对答案的贡献就是 $\frac 1 n\sum\limits_{i=1}^n x_i\cdot p_i$。因此 $\sum\limits_{i=1}^n x_i\cdot p_i$ 就是答案。

考虑求出 $x$。令 $f_i$ 表示 $i$ 的父结点,$S(i)$ 表示 $i$ 的祖先集合(包括自己),$\mathrm{size}(i)$ 表示 $i$ 的子树大小。

对每个结点 $v$,我们均能列出一个方程:

其中 $v=1$ 时的方程没有用。这样能得到 $n-1$ 个方程,有任意一组解。

再加上 $\sum\limits_{i=1}^n x_i=1$ 就有 $n$ 个方程了,这样就有唯一解了。

直接高斯消元解方程是 $O(n^3)$ 的。考虑到这是一个树上的问题,不难想到对每个 $i$ 找到系数 $k,b$ 使得 $x_i=kx_{fa_i}+b$ 然后向上传递。本题中,除了最后一个方程外,其它方程的常数项为 $0$,因此 $b=0$。

但是这 $n-1$ 个方程不仅和儿子结点有关,还和祖先结点有关,难以直接求解。那么我们考虑将和祖先结点有关的部分去掉。

令 $X_v=\sum_{i\in S(v)} \frac{1-p_i}{\mathrm{size}(i)}\cdot x_i$,则 $x_v=(X_v-X_{fa_v})\cdot\frac{\mathrm{size}(v)}{1-p_v}$。特别地,$X_1=\frac{x_1}{\mathrm{size}(1)}$。

将其代入原来的方程中,替换掉原来的 $x$,可以得到:

这个 $X$ 的方程组就只和每个结点的儿子结点有关了。我们对每个 $v$,从下到上解得满足 $X_v=kX_{f_v}$ 的 $k$。从而可以得到 $X_i$ 和 $X_1$ 的关系,再解得 $x_i$ 和 $x_1$ 的关系,最后通过 $x$ 的和为 $1$ 的条件即可解出所有的 $x$。

这可以在 $O(n)$ 的时间复杂度内完成。


接下来考虑加上 $x$ 的情况。

由于 $a$ 是一个排列,所以加上 $x$ 取模后仍是一个排列,那么必然存在一个 $i$ 满足 $a_i’ = 0$,且每个结点刚好满足一次。

由于轮数足够多,每只猴子一定会以 1 的概率在有限步内跳到满足 $a_i’=0$ 的结点的子树内,然后就不可能再上去了。因此只需要对这棵子树求解。

这个算法的复杂度为 $O(\sum\limits_{i=1}^n\mathrm{size(i)})$。由于树是随机父亲的方式生成的,因此它的复杂度是 $O(n\log n)$。

Code:

#include<cstdio>
#include<iostream>
using namespace std;
const int N=5e5+5;
int nxt[N],pr[N],n,sz[N],dfn[N],idx,a[N],ssz[N],A[N],fa[N],head[N];
double ans,p[N],q[N],x[N],X[N],K[N];
void dfs(int u){
    sz[u]=1,dfn[u]=++idx;
    for(int i=head[u];i;i=nxt[i])dfs(i),sz[u]+=sz[i];
}
void work(int r){
    const int L=dfn[r],R=dfn[r]+sz[r]-1;
    for(int i=R;i>=L;--i){
        p[i]=1.*((a[i]-A[r]+n)%n)/n;
        q[i]=ssz[i]/(1-p[i]);
        K[i]=q[i]-1;
    }
    for(int i=R;i>L;--i)K[i]=q[i]/K[i],K[fa[i]]-=(K[i]-1)*q[i]*p[i];
    X[L]=x[L]=1;
    for(int i=L+1;i<=R;++i)X[i]=X[fa[i]]*K[i],x[i]=X[i]-X[fa[i]];
    double s=0;
    for(int i=L;i<=R;++i)s+=(x[i]=x[i]*q[i]);
    x[1]=1/s;
    double res=0;
    for(int i=L;i<=R;++i)res+=x[i]*x[1]*p[i];
    if(ans<res)ans=res;
}
int main(){
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    cin>>n>>fa[1];
    for(int i=2;i<=n;++i)cin>>pr[i],nxt[i]=head[pr[i]],head[pr[i]]=i;
    dfs(1);
    for(int i=1;i<=n;++i)cin>>A[i],a[dfn[i]]=A[i],ssz[dfn[i]]=sz[i],fa[dfn[i]]=dfn[pr[i]];
    for(int i=1;i<=n;++i)
    work(i);
    printf("%.12f\n",ans);
    return 0;
}