介于网上唯一一篇 PGF 实在过于简略，且他的码风使我大受震撼，因此我决定自己写一篇
【大意】
初始给定序列 ({a_{n-1}}) 与 ({c_n}) ，从 (1) 到 (n) 构建一颗树：
对于当前点 (i) ，它有 (displaystyle {a_jover 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)，则显然 (kleq u)
- 当 (k=u) 时，不妨假设中间仍连着 ({w1, w2, cdots , wm}) 若干个点
则该情况的概率为 (displaystyle {a_uover s_{w1-1}}cdot {a_{w1}over s_{w2-1}}cdots {a_{wm}over s_{v-1}}={a_uover 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+2sum_{i=1}^m c_{wi})
考虑此时每个 (u<j<v) 都可以选或不选，故生成函数为 (displaystyle {a_ucdot x^{c_u+c_v}over s_{v-1}}prod_{j=u+1}^{v-1}(1+{a_jover 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_kover s_{w1-1}}cdot {a_{w1}over s_{w2-1}}cdots {a_{wm}over s_{u-1}})cdot ({a_kover s_{t1-1}}cdot {a_{t1}over s_{t2-1}}cdots {a_{tm}over s_{v-1}})={a_k^2over s_{u-1}cdot s_{v-1}}prod_{rin L_1cup L_2-{u, v}} {a_rover 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+2sum_{rin L_1cup L_2-{u, v}} c_r)
故同理考虑生成函数，理应为 (displaystyle {a_k^2cdot 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_jover s_{j-1}}x^{2c_j})cdot prod_{j=u+1}^{v-1}(1+{a_jover s_{j-1}}x^{2c_j}))
但这样的做法是有问题的，考虑一个链上小于 (u) 的数字 (q)，其他数字不变
当 (q) 出现在 (u) 所在的链上时，对概率会产生形如 (displaystyle {a_kover s_{q-1}}cdot {a_qover s_{u-1}}cdot {a_kover s_{v-1}}) 乘上其他数的贡献
而当 (q) 出现在 (v) 所在的链上时，路径长度不变，且对概率会产生形如 (displaystyle {a_kover s_{q-1}}cdot {a_qover s_{v-1}}cdot {a_kover s_{u-1}}) 乘上其他数的贡献
两者贡献相同，故应该将贡献由 (displaystyle (1+{a_qover s_{q-1}}x^{2c_q})) 改为 (displaystyle (1+{2a_qover s_{q-1}}x^{2c_q}))
故其生成函数应改为 (displaystyle {a_k^2cdot 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_jover s_{j-1}}x^{2c_j})cdot prod_{j=u+1}^{v-1}(1+{a_jover s_{j-1}}x^{2c_j}))
因此 (displaystyle P(x)=sum_{k=1}^{u-1}{a_k^2cdot 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_jover s_{j-1}}x^{2c_j})cdot prod_{j=u+1}^{v-1}(1+{a_jover s_{j-1}}x^{2c_j})+{a_ucdot x^{c_u+c_v}over s_{v-1}}prod_{j=u+1}^{v-1}(1+{a_jover s_{j-1}}x^{2c_j}))
这个式子太复杂了，我们记 (displaystyle p1[n]=sum_{i=1}^{n-1}(1+{a_jover s_{j-1}}x^{2c_j}), p2[n]=sum_{i=1}^{n-1}(1+{2a_jover s_{j-1}}x^{2c_j}))
代入可化简得
(egin{aligned} P(x)&=sum_{k=1}^{u-1}{a_k^2cdot 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_jover s_{j-1}}x^{2c_j})cdot prod_{j=u+1}^{v-1}(1+{a_jover s_{j-1}}x^{2c_j})+{a_ucdot x^{c_u+c_v}over s_{v-1}}prod_{j=u+1}^{v-1}(1+{a_jover 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_nover 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) ight>) 同步维护它和它的导函数
然后定义二元运算 (left<f(x), f'(x) ight> imes left<g(x), g'(x) ight>=left<f(x)g(x), f'(x)g(x)+f(x)g'(x) ight>) ，那么我们就直接算出了它们乘积的导函数
很快我们发现，这个时候直接维护 (left<f(1), f'(1) ight>) ，也可以计算后直接得出乘积的单点函数值，以及乘积导函数的单点函数值
于是我们过程中遇到的所有函数全部维护 (left<f(1), f'(1) ight>) ，然后再按导数运算性质重新定义运算
最后算出的答案即是 (left<P(1), P'(1) ight>)
对于所有询问 (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)<<"
";
cout.flush();
return 0;
}
【附注】
这里还涉及到一个细节，即是无法保证 (P(1)=1) ，但本人暂时无法证明，只能从列式角度感觉它是“完美”的，因此保证计算结果为 (1)
于是本人在代码上加上了 assert( (res.fi+P)%P==1 ) ，实践证明是没有问题的
等后续有办法证明时，会回来证明