1. 等价转化

考虑一个等价的游走方式：

1. 如果所在的结点非叶子，以 $\frac1{\deg_x}$ 的概率走到 $x$ 的每个子结点。
2. 如果所在的结点是叶子，那么移到它的父亲，接下来以 $\frac{1}{dep_x}$ 的概率留在 $x$，否则继续移到 $fa_x$，以此类推。注意连续向上跳的过程不计入步数。

2. 求解 $out_x$

设 $x\ne 1$，$out_x$ 表示从 $x$ 出发，第一次移出 $x$ 的子树的期望步数，$A_x$ 表示 $x$ 的祖先结点构成的集合。若 $x$ 为叶子，则 $out_x=1$，否则 $out_x=\frac{dep_x}{dep_x-1}\left(1+\frac{1}{\deg_x}\sum_{v\in Son_x}out_v\right)$。

自下而上求出每个 $out_x$。

3. 求解 $E_x,down_x$

设 $E_x$ 表示从根出发，第一次经过 $x$ 的期望时间。$down_x$ 表示从 $fa_x$ 出发第一次到 $x$ 的期望时间，设 $f=fa_x$。

则有 $E_x=E_{f}+down_x$。

$down_x=\frac1{\deg_f}\left(1+\sum_{v\in Son_f,v\ne x}\left(1+out_v+\frac{1}{dep_f}\sum_{y\in A_f}\left(E_f-E_y\right)+down_x\right)\right)$。

自上而下联立去解即可。

4. 求解答案

设 $l=LCA(x,y)$。

若 $x=l$，答案为 $E_y-E_x$。

否则，考虑 $x$ 到 $y$ 的路径可以拆为三段：

1. 从 $x$ 到 $l$，设这一部分的期望时间为 $g_{fa_x,l}+out_x$。
2. 从 $l$ 继续走到 $l$ 的某个祖先 $l'$。期望时间为 $0$。
3. 从 $l'$ 到 $y$。期望时间为 $E_y-E_{l'}$。

设 $g_{s,l}$ 表示从 $s$ 出发向上走，第一次经过 $l$ 的期望时间。则 $g_{s,s}=0$，$g_{s,l}=\frac1{dep_s}out_s+g_{fa_s,l}$。展开，得到 $g_{s,l}=\sum_{u\in Path(s,l),u\ne l}\frac{out_u}{dep_u}$。

将三部分的期望时间相加就是答案。

5. code

#include<bits/stdc++.h>
bool Mbg;
using namespace std;
#define vec vector
#define pb push_back
#define eb emplace_back
#define int long long
const int mod=1e9+7;
void Add(int &x,const int &y){x=x+y<mod?x+y:x+y-mod;}
void Dec(int &x,const int &y){x=x>=y?x-y:x-y+mod;}
int fpm(int x,int y){
    int ans=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)ans=1ll*x*ans%mod;return ans;
}

int n,q,iv[200010];
vec<int> Son[200010];
int out[200010],dep[200010],down[200010],E[200010],deg[200010];
int g[200010],fa[200010];
int anc[500010][20];
void dfs1(int u){
    anc[u][0]=fa[u];
    for(int i=1;i<20;i++)anc[u][i]=anc[anc[u][i-1]][i-1];
    dep[u]=dep[fa[u]]+1;
    for(auto v:Son[u]){
        ++deg[u];
        dfs1(v);
    }
    int x=u;
    if(x!=1){
        if(!deg[x]){
            out[x]=1;
        }else{
            for(auto v:Son[u]){
                (out[x]+=out[v])%=mod;
            }
            out[x]=(out[x]*iv[deg[x]]+1)%mod*dep[x]%mod*iv[dep[x]-1]%mod;
        }
    }
}
int s[200010],sE[200010];
void dfs2(int u){
    int f=fa[u];
    g[u]=(g[fa[u]]+out[u]*iv[dep[u]])%mod;
    int x=u;

    down[x]=(s[f]-out[u]+mod)*iv[deg[f]]+1+(deg[f]-1)*iv[deg[f]]%mod*(E[f]-iv[dep[f]]*sE[f]%mod+mod)%mod;
    (down[x])%=mod;
    (down[x]*=deg[f])%=mod;

    E[x]=(E[f]+down[x])%mod;
    sE[x]=(sE[f]+E[x])%mod;
    for(auto v:Son[u]){
        (s[u]+=out[v])%=mod;
    }
    for(auto v:Son[u]){
        dfs2(v);
    }
}

int LCA(int u,int v){
    if(dep[u]<dep[v])swap(u,v);
    for(int i=19;~i;i--){
        if((dep[u]-dep[v])>>i&1)u=anc[u][i];
    }
    if(u==v)return u;
    for(int i=19;~i;i--){
        if(anc[u][i]!=anc[v][i])u=anc[u][i],v=anc[v][i];
    }
    return anc[u][0];
}

void work(){
    cin>>n>>q;
    for(int i=2;i<=n;i++){
        cin>>fa[i];
        Son[fa[i]].eb(i);
    }
    for(int i=1;i<=n;i++)iv[i]=fpm(i,mod-2);
    dfs1(1);
    dfs2(1);
    while(q--){
        int x,y;
        cin>>x>>y;
        int l=LCA(x,y);
        if(x==y)cout<<"0\n";
        else if(x==l){
            cout<<(E[y]-E[x]+mod)%mod<<'\n';
        }else{
            int ans=(out[x]+g[fa[x]]-g[l]+mod)%mod;
            (ans+=E[y]-sE[l]*iv[dep[l]]%mod+mod)%=mod;
            cout<<ans<<'\n';
        }
    }
}
bool Med;
signed main(){
    ios::sync_with_stdio(0),
    cin.tie(0),cout.tie(0);
    int T=1;while(T--)work();
    // cerr<<"Time: "<<clock()<<" ms;\n";
    // cerr<<"Memory: "<<abs(&Med-&Mbg)/1024.0/1024.0<<" MiB.\n";
}

(Analysis by Arvith Vijayram)

Let $up_x$ be the expected time for Bessie to reach an ancestor of $x$ for the first time when starting at $x$. If $x$ is a leaf node, then $up_x = 1$. Otherwise, the expected time to either reach $x$ again or an ancestor of $x$ is the average of $up_c$ over all children $c$ of $x$ plus 1. Since all ancestors of $x$ are equally likely to be the first one visited in this scenario, the probability we land on $x$ again is $\frac{1}{d_x+1}$, where $d_x$ is the number of ancestors of $x$. Thus, we have

$$up_x = 1 + \frac{1}{|C_x|} \sum_{c \in C_x} up_c + \frac{1}{d_x+1} up_x$$

Here, $C_x$ is the set of children of node $x$. We can solve for $up_x$, giving us

$$up_x = \left(1 + \frac{1}{d_x}\right) \left(1 + \frac{1}{|C_x|} \sum_{c \in C_x} up_c\right)$$

We can compute $up_x$ for all nodes $x$ (except the root) with a DFS in $O(N)$.

Subtask 1:

Let $root_x$ represent the expected time to reach node 1 starting at node $x$. We have $root_1 = 0$. Additionally, letting $A_x$ denote the set of ancestors of $x$, we have

$$root_x = up_x + \frac{1}{d_x} \sum_{a \in A_x} root_a$$

We can compute $root_x$ for all $x$ in $O(N)$ with a DFS from node 1 by propagating the sum $\sum_{a \in A_x} root_a$ down the tree.

Subtask 2:

Let $E_y$ be the expected time to reach node $y$ from node 1. Define $down_y$ be the expected time to reach node $y$ from $par_y$, the parent of $y$, and let $down_1 = 0$. We have $E_y = \sum_{a \in A_y \cup \{y\}} down_a$.

If we start from $par_y$, we either take the correct step to $y$ with chance $\frac{1}{|C_{par_y}|}$ but otherwise have to return to $par_y$ or an ancestor of $par_y$. Thus,

$$down_y = 1 + \frac{|C_{par_y}| - 1}{|C_{par_y}|} \left( \underset{c \in C_{par_y} \setminus \{y\}}{\text{avg}} up_c + \underset{a \in A_y}{\text{avg}} (E_{par_y} - E_a) + down_y \right)$$

$E_{par_y} - E_a$ represents the time it takes to reach $par_y$ from $a$ since every path from 1 to $par_y$ goes through $a$. We can rearrange the equation to the following form.

$$down_y = |C_{par_y}| + \sum_{c \in C_{par_y} \setminus \{y\}} up_c + (|C_{par_y}| - 1) \left( E_{par_y} - \frac{1}{d_y} \sum_{a \in A_y} E_a \right)$$

Notice that $down_y$ only depends on the values $E_a$ for all ancestors $a$ of $y$. By maintaining $E_{par_y}$ and $\sum_{a \in A_y} E_a$, we can compute $down_y$ and $E_y$ for all nodes $y$ with a DFS in $O(N)$.

Answering Queries:

Given nodes $x, y$, we want to determine the expected time to reach $y$ starting at $x$. Let $S$ be the set of common ancestors of $x$ and $y$ and let $\ell \in S$ be the lowest common ancestor (LCA). Any path from $x$ to $y$ must go through a node in $S$. Let $w$ be the first node in $S$ we visit when starting at node $x$. We can split the path into $x \to w$ and $w \to y$. By linearity of expectation, the answer is the average expected time from $x \to w$ over all $w \in S$ plus the average expected time from $w \to y$ over all $w \in S$. The latter is $\frac{1}{|S|} \sum_{w \in S} (E_y - E_w) = E_y - \frac{1}{|S|} \sum_{w \in S} E_w$ by linearity of expectation since the path from 1 to $y$ must pass through $w$.

Now, we want to compute the average expected time from $x \to w$ over all $w \in S$. Since all $w \in S$ are equally likely, we need to find the expected time to reach any $w \in S$. Let $\ell = x_0, x_1, x_2, \dots, x_m = x$ be the nodes on the path from $\ell$ to $x$. Define $reach_i$ as the expected time to reach a node in $S$ when starting at $x_i$. Our goal is to compute $reach_m$.

We can observe the following relation:

$$reach_i = up_{x_i} + \frac{1}{d_{x_i}} \sum_{j=1}^{i-1} reach_j$$

This is because it takes expected time $up_{x_i}$ to move to an ancestor of $x_i$, and we have an $\frac{|S|}{d_x}$ chance to land in $S$ and an $\frac{1}{d_x}$ chance each to land on $x_1, x_2, \dots, x_{i-1}$.

Subtask 3:

In this subtask, we use the fact that the expected maximum depth of the tree is $O(\log N)$. Thus, we can compute $reach_m$ in $O(m^2)$ or $O(m)$ per query for a total time complexity of $O(N + Q \log^2 N)$ or $O(N + Q \log N)$, respectively.

Full solution:

We claim that

$$reach_i = up_{x_i} + \sum_{j=1}^{i-1} \frac{up_{x_j}}{d_{x_j} + 1}$$

This is because for each $x_j$, the probability that it is visited before any of its ancestors is $\frac{1}{d_{x_j}+1}$.

Thus, to compute $reach_m$ we just need to sum $\frac{up_z}{d_z+1}$ for all $x_j$ on the path from $\ell$ to $y$, both exclusive. For each $v$, we can precompute the sum of $\frac{up_z}{d_z+1}$ over all nodes $z$ on the path from node 1 to node $v$ in $O(N)$. This allows us to compute the sum along a path by taking differences of two of these sums. We now have a solution in $O(N \log N + Q \log N)$, though this may vary slightly depending on implementation of LCA.

#include <bits/stdc++.h>
using namespace std;

using ll = long long;
const int maxn = 2e5+5;
const ll M = 1e9+7;

ll modpow(ll x, ll p) {
    ll a = 1;
    while (p) {
        if (p & 1) a = a * x % M;
        x = x * x % M;
        p /= 2;
    }
    return a;
}

ll inv(ll x) {
    assert(x != 0);
    return modpow(x, M-2);
}

ll par[maxn];
vector<ll> childs[maxn];

ll d[maxn], up[maxn], e[maxn], esum[maxn], upd[maxn];

ll dfs1(int x, ll dep) {
    d[x] = dep;
    ll cx = childs[x].size();
    if (cx == 0) {
        up[x] = 1;
        return up[x];
    }

    ll upsum = 0;
    for (auto c : childs[x]) {
        upsum += dfs1(c, dep+1);
    }
    upsum %= M;

    if (dep > 0) {
        up[x] = (1 + inv(dep)) * (1 + inv(cx) * upsum % M) % M;
    }
    return up[x];
}

void dfs2(int x, ll cp, ll upcsum, ll ep, ll epsun, ll updsum) {
    ll downx = 0;
    if (x != 0) {
        downx = (cp + upcsum + (cp - 1) * (ep - inv(d[x]) * epsun % M)) % M;
        if (downx < 0) downx += M;
    }
    e[x] = (ep + downx) % M;
    esum[x] = (epsun + e[x]) % M;
    if (x > 0) updsum = (updsum + up[x] * inv(d[x]+1)) % M;
    upd[x] = updsum;

    ll cx = childs[x].size();
    ll upsum = 0;
    for (auto c : childs[x]) {
        upsum += up[c];
    }

    for (auto c : childs[x]) {
        dfs2(c, cx, (upsum + M - up[c]) % M, e[x], esum[x], upd[x]);
    }
}

int main() {
    int n, q; cin >> n >> q;
    for (int i = 1; i < n; i++) {
        int p; cin >> p; p--;
        childs[p].push_back(i);
        par[i] = p;
    }

    dfs1(0, 0);
    dfs2(0, 0, 0, 0, 0, 0);

    vector<vector<int>> jmp(18, vector<int>(n, -1));
    for (int i = 1; i < n; i++) jmp[0][i] = par[i];
    for (int j = 1; j < 18; j++) {
        for (int i = 0; i < n; i++) {
            int nx = jmp[j-1][i];
            if (nx != -1) nx = jmp[j-1][nx];
            jmp[j][i] = nx;
        }
    }

    auto jump = [&](int x, int l) -> int {
        for (int j = 17; j >= 0; j--) {
            if (x == -1) return x;
            if (l & (1 << j)) x = jmp[j][x];
        }
        return x;
    };

    while (q--) {
        int x, y; cin >> x >> y; x--; y--;

        int l;
        int nx = x, ny = y;
        if (d[nx] > d[ny]) nx = jump(nx, d[nx]-d[ny]);
        if (d[ny] > d[nx]) ny = jump(ny, d[ny]-d[nx]);
        if (nx == ny) l = nx;
        else {
            for (int j = 17; j >= 0; j--) {
                if (jmp[j][nx] != jmp[j][ny]) {
                    nx = jmp[j][nx];
                    ny = jmp[j][ny];
                }
            }
            assert(par[nx] == par[ny]);
            l = par[nx];
        }

        if (l == x) {
            cout << (e[y]-e[x]+M)%M << endl;
        } else {
            ll reachm = (up[x] + upd[par[x]] - upd[l] + M) % M;
            ll wtoy = (e[y] - (inv(d[l]+1) * esum[l] % M) + M) % M;
            cout << (reachm + wtoy) % M << endl;
        }
    }
}