题解 cf gym 102979 E Expected Distance
介于网上唯一一篇 PGF 实在过于简略，且他的码风使我大受震撼，因此我决定自己写一篇
【大意】
初始给定序列 \(\{a_{n-1}\}\) 与 \(\{c_n\}\) ，从 \(1\) 到 \(n\) 构建一颗树：
对于当前点 \(i\) ，它有 \(\displaystyle {a_j\over \sum_{t=1}^{i-1} a_t}\) 的概率连向 \(j\) ，当该事件发生时， \((i, j)\) 的边长为 \((c_i+c_j)\)
现在 \(Q\) 次询问，每次询问输出 \(u\) 到 \(v\) 的距离期望
【分析】
先定义概率型生成函数 PGF ：
假定对于离散型随机事件 \(X\) ，有 \(P(X=i)=p_i\) ，则其 PGF 定义为 \(\displaystyle P(x)=\sum_{i=0}^\infty P(X=i)x^i=\sum_{i=0}^\infty p_ix^i\)
根据概率的性质，显然有 \(P(1)=1\) 时，\(\displaystyle \sum_{i=0}^\infty p_i=1\) ，则 \(\displaystyle E(X)=\sum_{i=0}^\infty ip_i=\sum_{i=1}^\infty ip_i=P'(1)\)
对于该题，由于询问内容是距离的期望，因此我们假设 \((u, v)\) 距离的 PGF 为 \(\displaystyle P(x)=\sum_{i=0}^\infty P[dis(u, v)=i]x^i\) ，则答案为求解 \(P'(1)\)
不妨设 \(u<v\) ，一边相连的两点中，数值较大的为孩子节点，且为了方便，记 \(\displaystyle s_n=\sum_{i=1}^n a_i\)
考虑枚举 \((u, v)\) 的最近公共祖先 \(k\)，则显然 \(k\leq u\)
- 当 \(k=u\) 时，不妨假设中间仍连着 \(\{w1, w2, \cdots , wm\}\) 若干个点
则该情况的概率为 \(\displaystyle {a_u\over s_{w1-1}}\cdot {a_{w1}\over s_{w2-1}}\cdots {a_{wm}\over s_{v-1}}={a_u\over s_{v-1}}\prod_{i=1}^m {a_{wi}\over s_{wi-1}}\) ，此时的路径长为 \(\displaystyle c_u+2c_{w1}+2c_{w2}+\cdots +2c_{wm}+c_v=c_u+c_v+2\sum_{i=1}^m c_{wi}\)
考虑此时每个 \(u<j<v\) 都可以选或不选，故生成函数为 \(\displaystyle {a_u\cdot x^{c_u+c_v}\over s_{v-1}}\prod_{j=u+1}^{v-1}(1+{a_j\over s_{j-1}}x^{2c_j})\)
- 当 \(k<u\) 时，不妨假设两条链分别为 \(L_1=\{k, w1, w2, \cdots , wm, u\}, L_2=\{k, t1, t2, \cdots, tn, v\}\)
则该情况的概率为 \(\displaystyle ({a_k\over s_{w1-1}}\cdot {a_{w1}\over s_{w2-1}}\cdots {a_{wm}\over s_{u-1}})\cdot ({a_k\over s_{t1-1}}\cdot {a_{t1}\over s_{t2-1}}\cdots {a_{tm}\over s_{v-1}})={a_k^2\over s_{u-1}\cdot s_{v-1}}\prod_{r\in L_1\cup L_2-\{u, v\}} {a_r\over s_{r-1}}\) ，此时的路径长为 \(\displaystyle (c_k+2c_{w1}+2c_{w2}+\cdots +2c_{wm}+c_u)+(c_k+2c_{t1}+2c_{t2}+\cdots +2c_{tn}+c_v)=c_u+c_v+2\sum_{r\in L_1\cup L_2-\{u, v\}} c_r\)
故同理考虑生成函数，理应为 \(\displaystyle {a_k^2\cdot x^{2c_k}\cdot x^{c_u+c_v}\over s_{u-1}\cdot s_{v-1}}\prod_{j=k+1}^{u-1}(1+{a_j\over s_{j-1}}x^{2c_j})\cdot \prod_{j=u+1}^{v-1}(1+{a_j\over s_{j-1}}x^{2c_j})\)
但这样的做法是有问题的，考虑一个链上小于 \(u\) 的数字 \(q\)，其他数字不变
当 \(q\) 出现在 \(u\) 所在的链上时，对概率会产生形如 \(\displaystyle {a_k\over s_{q-1}}\cdot {a_q\over s_{u-1}}\cdot {a_k\over s_{v-1}}\) 乘上其他数的贡献
而当 \(q\) 出现在 \(v\) 所在的链上时，路径长度不变，且对概率会产生形如 \(\displaystyle {a_k\over s_{q-1}}\cdot {a_q\over s_{v-1}}\cdot {a_k\over s_{u-1}}\) 乘上其他数的贡献
两者贡献相同，故应该将贡献由 \(\displaystyle (1+{a_q\over s_{q-1}}x^{2c_q})\) 改为 \(\displaystyle (1+{2a_q\over s_{q-1}}x^{2c_q})\)
故其生成函数应改为 \(\displaystyle {a_k^2\cdot x^{2c_k}\cdot x^{c_u+c_v}\over s_{u-1}\cdot s_{v-1}}\prod_{j=k+1}^{u-1}(1+{2a_j\over s_{j-1}}x^{2c_j})\cdot \prod_{j=u+1}^{v-1}(1+{a_j\over s_{j-1}}x^{2c_j})\)
因此 \(\displaystyle P(x)=\sum_{k=1}^{u-1}{a_k^2\cdot x^{2c_k}\cdot x^{c_u+c_v}\over s_{u-1}\cdot s_{v-1}}\prod_{j=k+1}^{u-1}(1+{2a_j\over s_{j-1}}x^{2c_j})\cdot \prod_{j=u+1}^{v-1}(1+{a_j\over s_{j-1}}x^{2c_j})+{a_u\cdot x^{c_u+c_v}\over s_{v-1}}\prod_{j=u+1}^{v-1}(1+{a_j\over s_{j-1}}x^{2c_j})\)
这个式子太复杂了，我们记 \(\displaystyle p1[n]=\sum_{i=1}^{n-1}(1+{a_j\over s_{j-1}}x^{2c_j}), p2[n]=\sum_{i=1}^{n-1}(1+{2a_j\over s_{j-1}}x^{2c_j})\)
代入可化简得
\(\begin{aligned} P(x)&=\sum_{k=1}^{u-1}{a_k^2\cdot x^{2c_k}\cdot x^{c_u+c_v}\over s_{u-1}\cdot s_{v-1}}\prod_{j=k+1}^{u-1}(1+{2a_j\over s_{j-1}}x^{2c_j})\cdot \prod_{j=u+1}^{v-1}(1+{a_j\over s_{j-1}}x^{2c_j})+{a_u\cdot x^{c_u+c_v}\over s_{v-1}}\prod_{j=u+1}^{v-1}(1+{a_j\over s_{j-1}}x^{2c_j}) \\\\&=(\sum_{k=1}^{u-1}{a_k^2x^{2c_k}\over s_{u-1}s_{v-1}}\cdot {p_2[u]\over p_2[k+1]}\cdot {p1[v]\over p1[u+1]}+{a_n\over s_{v-1}}\cdot {p_1[v]\over p_1[u+1]})\cdot x^{c_u+c_v} \\\\&={p_1[v]\cdot x^{c_u+c_v}\over s_{v-1}p_1[u+1]}({p2[u]\over s_{u-1}}\sum_{k=1}^{u-1}{a_k^2x^{2c_k}\over p_2[k+1]}+a_n) \end{aligned}\)
于是我们除了预处理 \(p1[n], p2[n]\) ，在额外预处理 \(\displaystyle sum[n]=\sum_{i=1}^{n-1}{a_i^2x^{2c_i}\over p_2[i+1]}\)
则可通过 \(\displaystyle P(x)={p_1[v]\cdot x^{c_u+c_v}\over s_{v-1}p_1[u+1]}({p2[u]\over s_{u-1}}\cdot sum[u]+a_n)\) 求解答案了
但是上述式子有的是多项式，有的是常数，最后计算 \(P(x)\) 再求导算出 \(P'(1)\) 过于复杂了
这边考虑一个性质：假设对于一个多项式 \(f(x)\) ，我们用二元组 \(\left<f(x), f'(x)\right>\) 同步维护它和它的导函数
然后定义二元运算 \(\left<f(x), f'(x)\right>\times \left<g(x), g'(x)\right>=\left<f(x)g(x), f'(x)g(x)+f(x)g'(x)\right>\) ，那么我们就直接算出了它们乘积的导函数
很快我们发现，这个时候直接维护 \(\left<f(1), f'(1)\right>\) ，也可以计算后直接得出乘积的单点函数值，以及乘积导函数的单点函数值
于是我们过程中遇到的所有函数全部维护 \(\left<f(1), f'(1)\right>\) ，然后再按导数运算性质重新定义运算
最后算出的答案即是 \(\left<P(1), P'(1)\right>\)
对于所有询问 \(O(1)\) 直接输出结果即可
【代码】
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
typedef double db;
#define fi first
#define se second
const int MAXN=3e5+10, P=1e9+7;
inline ll fpow(ll a,ll x) { ll ans=1; for(;x;x>>=1,a=a*a%P) if(x&1) ans=ans*a%P; return ans; }
inline ll inv(ll a) { return fpow(a, P-2); }
inline pii operator + (const pii &a, const pii &b) { return pii( (a.fi+b.fi)%P, (a.se+b.se)%P); }
inline pii operator - (const pii &a, const pii &b) { return pii( (a.fi-b.fi)%P, (a.se-b.se)%P); }
inline pii operator * (const pii &a, const pii &b) { return pii( (ll)a.fi*b.fi%P, ((ll)a.fi*b.se+(ll)a.se*b.fi)%P ); }
inline pii operator / (const pii &a, const pii &b) { int x=inv(b.fi); return pii( (ll)a.fi*x%P , ((ll)a.se*b.fi-(ll)a.fi*b.se)%P*x%P*x%P ); }
inline pii operator + (const pii &a, int b) { return pii( (a.fi+b)%P, a.se ); }
inline pii operator - (const pii &a, int b) { return pii( (a.fi-b)%P, a.se ); }
inline pii operator * (const pii &a, int b) { return pii( (ll)a.fi*b%P, (ll)a.se*b%P ); }
inline pii operator / (const pii &a, int b) { int x=inv(b); return pii( (ll)a.fi*x%P, (ll)a.se*x%P ); }
int n, q, a[MAXN], c[MAXN], s[MAXN];
pii prod1[MAXN], prod2[MAXN], sumit[MAXN];
inline void init() {
cin>>n>>q;
for(int i=1; i<n; ++i) cin>>a[i], s[i]=(s[i-1]+a[i])%P; s[n]=s[n-1];
for(int i=1; i<=n; ++i) cin>>c[i];
prod1[1]=prod2[1]=prod1[2]=prod2[2]=pii(1, 0);
for(int i=3; i<=n; ++i) {
int x=a[i-1]*inv(s[i-2])%P;
pii tmp(1, (c[i-1]+c[i-1])%P);
prod1[i]=tmp*x+1;
prod2[i]=tmp*(x+x)+1;
prod1[i]=prod1[i-1]*prod1[i];
prod2[i]=prod2[i-1]*prod2[i];
}
sumit[1]=pii(0, 0);
for(int i=2; i<=n; ++i) {
int x=(ll)a[i-1]*a[i-1]%P;
sumit[i]=pii(1, (c[i-1]+c[i-1])%P)*x/prod2[i];
sumit[i]=sumit[i]+sumit[i-1];
}
}
inline int ans(int u, int v) {
if(u==v) return 0;
if(u>v) swap(u, v);
pii res=sumit[u]*prod2[u]/s[u-1]+a[u];
res=res*prod1[v]*pii(1, (c[u]+c[v])%P)/prod1[u+1]/s[v-1];
assert( (res.fi+P)%P==1 );
return (res.se+P)%P;
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
init();
int u, v;
while(q--&&cin>>u>>v)
cout<<ans(u, v)<<"\n";
cout.flush();
return 0;
}
【附注】
这里还涉及到一个细节，即是无法保证 \(P(1)=1\) ，但本人暂时无法证明，只能从列式角度感觉它是“完美”的，因此保证计算结果为 \(1\)
于是本人在代码上加上了 assert( (res.fi+P)%P==1 ) ，实践证明是没有问题的
等后续有办法证明时，会回来证明
