T3
statement
有一个有权完全图，每次等概率选择一个点对 \((u,v)\)，把这条边的权值 \(v=(v+1)\bmod 3\),\(q\) 次询问，求 \(t\) 次操作后的图的所有生成树的边权乘积的和的期望。
\(n\leq 50,q\leq 10^4,t\leq 10^{18}\)
solution
首先经典套路，把每条边看成一个二元多项式，求所有生成树的多项式乘积和 \(T\)，那么 \([x^iy^j]T\) 就是 \(i\) 条 \(1\) 边，\(j\) 条 \(2\) 边的生成树个数，这个可以用二元拉格朗日插值求出来。
可以 \(dp\) 出来一个 \(f_{i,j}\)表示 假如最后有 \(i\) 条边被 \(+1\)，\(j\) 条边被 \(+2\)，\(n-1-i-j\) 条边被 \(+0\)，剩下 \(\binom{n}{2}-n+1\) 条边无限制的所有生成树的权值和。
现在我们就要把 \(t\) 次操作分配给这些东西，考虑列出 \(EGF\) ，并使用单位根反演，设 \(3\) 次单位根为 \(w\)。
假设我们要求出
\[\begin{align} F_k(x) & = \sum\limits_{i=0}[i\bmod 3=k]\frac{x^i}{i!} \\ & = \sum\limits_{i=0}[3|(i-k)]\frac{x^i}{i!} \\ & = \frac{1}{3}\sum\limits_{i=0}\sum\limits_{j=0}^2w^{j(i-k)}\frac{x^i}{i!} \\ & = \frac{1}{3}\sum\limits_{j=0}^2w^{-jk}\sum\limits_{i=0}\frac{(w^j x)^i}{i!} \\ & = \frac{1}{3}\sum\limits_{j=0}^2w^{-jk}e^{w^j x} \\ \end{align} \]-
边权不变
\[F_0=\frac{1}{3} (e ^ {w^2 x} + e ^ {w x} + e ^x) \] -
边权加一
\[F_1=\frac{1}{3}(w^2 e^{wx}+we^{w^2x}+e^x) \] -
边权加二
\[F_2=\frac{1}{3}(we^{wx}+w^2e^{w^2x}+e^x) \]
那么我们要求
\[[x^t]F_0^{n-1-i-j}F_1^{i}F_2^{j}(e^x)^{\binom{n}{2}-n+1} \]因为 \(w^2=-w-1\)，所以上述 \(EGF\) 一定能表示成 \(\sum_{a,b}[x^t]e^{(aw+b)x}g_{a,b}\) 的形式，\((a,b)\) 只有 \(O(n^2)\) 个，所以可以 \(O(n^2\log t)\) 计算，可以用光速幂去掉 $\log $，总复杂度 \(O(n^5+qn^2)\)。
code
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 55;
const int mod = 1e9+9;
const int i3=(mod*2+1)/3;
inline void myplus(int &x,int y){x+=y;if(x>=mod) x-=mod;}
inline void reduce(int &x,int y){x-=y;if(x<0) x+=mod;}
LL Pow(LL a,LL b)
{
LL res=1;
while(b)
{
if(b&1)res=res*a%mod;
a=a*a%mod;
b>>=1;
}
return res;
}
inline LL Inv(LL x){return Pow(x,mod-2);}
const LL w=Pow(13,(mod-1)/3),w2=w*w%mod;
int n,m,q;
int e[N][N];
int inv[N],fac[N],ifac[N];
void pre()
{
inv[1]=fac[1]=fac[0]=ifac[0]=ifac[1]=1;
for(int i=2;i<=n;i++)
{
inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
fac[i]=1ll*fac[i-1]*i%mod;
ifac[i]=1ll*ifac[i-1]*inv[i]%mod;
}
}
int p[N][N];
namespace det
{
int d[N][N];
int f[N][N];
int det(int x)
{
int res=1;
for(int i=1;i<=x;i++)
{
if(!d[i][i])
{
int id=i;
for(int j=i+1;j<=x;j++) if(d[j][i]) id=j;
swap(d[id],d[i]);res=-res;
}
LL inv=Inv(d[i][i]);
for(int j=i+1;j<=x;j++)
{
int t=1ll*d[j][i]*inv%mod;
for(int k=i;k<=x;k++) reduce(d[j][k],1ll*d[i][k]*t%mod);
}
}
res+=mod;
for(int i=1;i<=x;i++) res=1ll*res*d[i][i]%mod;
return res;
}
int B[N][N];
int F[N],G[N];
void init()
{
F[0]=1;
for(int i=1;i<=n;i++)
{
for(int j=i;j>=0;j--)
{
F[j]=1ll*F[j]*(mod-i)%mod;
if(j) myplus(F[j],F[j-1]);
}
}
}
void preinit(int x)
{
for(int i=1;i<=n;i++)
{
int w=f[x][i];
for(int j=1;j<=n;j++) if(j!=i) w=1ll*w*(i<j?mod-inv[j-i]:inv[i-j])%mod;
for(int j=n-1;j>=0;j--) G[j]=(F[j+1]+1ll*G[j+1]*i)%mod;
for(int j=0;j<n;j++) myplus(B[x][j],1ll*w*G[j]%mod);
}
}
void lagrange()
{
for(int i=1;i<=n;i++)
{
int w=1;
for(int j=1;j<=n;j++) if(j!=i) w=1ll*w*(i<j?mod-inv[j-i]:inv[i-j])%mod;
for(int j=n-1;j>=0;j--) G[j]=(F[j+1]+1ll*G[j+1]*i)%mod;
for(int j=0;j<n;j++) for(int k=0;k<n;k++)myplus(p[j][k],1ll*B[i][k]*G[j]%mod*w%mod);
}
}
void pre()
{
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
{
memset(d,0,sizeof(d));
for(int u=1;u<=n;u++)
for(int v=u+1;v<=n;v++)
{
int val=(e[u][v]==0?1:(e[u][v]==1?i:j));
myplus(d[u][u],val);myplus(d[v][v],val);
reduce(d[u][v],val);reduce(d[v][u],val);
}
f[i][j]=det(n-1);
}
init();
for(int i=1;i<=n;i++) preinit(i);
lagrange();
}
}
int f[N][N];
namespace F
{
int g[N][N];
void pre()
{
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
if(p[i][j])
{
memset(g,0,sizeof(g));
g[0][0]=p[i][j];
for(int k=1;k<=i;k++)for(int l=n-1;l>=0;l--)
if(l) myplus(g[l][0],1ll*g[l-1][0]*2%mod);
for(int k=1;k<=j;k++)for(int l=n-1;l>=0;l--)for(int o=n-1-l;o>=0;o--)
{
myplus(g[l][o],g[l][o]);
if(o) myplus(g[l][o],g[l][o-1]);
}
for(int k=1;k<=n-1-i-j;k++)for(int l=n-1;l>=0;l--)for(int o=n-1-l;o>=0;o--)
{
g[l][o]=0;
if(l) myplus(g[l][o],g[l-1][o]);
if(o) myplus(g[l][o],1ll*g[l][o-1]*2%mod);
}
for(int k=0;k<n;k++)for(int l=0;l+k<n;l++)myplus(f[k][l],g[k][l]);
}
}
}
int g[N*2][N*2];
namespace G
{
int f[N*2][N*2],tf[N*2][N*2];
void s1()
{
memset(tf,0,sizeof(tf));
for(int i=1;i<2*n;i++) for(int j=1;j<2*n;j++) if(f[i][j])
{
myplus(tf[i][j+1],w2*f[i][j]%mod);
myplus(tf[i-1][j-1],w*f[i][j]%mod);
myplus(tf[i+1][j],f[i][j]);
}
memcpy(f,tf,sizeof(f));
}
void s2()
{
memset(tf,0,sizeof(tf));
for(int i=1;i<2*n;i++) for(int j=1;j<2*n;j++) if(f[i][j])
{
myplus(tf[i][j+1],w*f[i][j]%mod);
myplus(tf[i-1][j-1],w2*f[i][j]%mod);
myplus(tf[i+1][j],f[i][j]);
}
memcpy(f,tf,sizeof(f));
}
void s0()
{
memset(tf,0,sizeof(tf));
for(int i=1;i<2*n;i++) for(int j=1;j<2*n;j++) if(f[i][j])
{
myplus(tf[i][j+1],f[i][j]);myplus(tf[i-1][j-1],f[i][j]);myplus(tf[i+1][j],f[i][j]);
}
memcpy(f,tf,sizeof(f));
}
void pre()
{
LL pw=Pow(i3,n-1);
for(int i=0;i<n;i++)for(int j=0;j<n;j++) if(::f[i][j])
{
memset(f,0,sizeof(f));
f[n][n]=pw*::f[i][j]%mod;
for(int k=1;k<=n-1-i-j;k++) s0();
for(int k=1;k<=i;k++) s1();
for(int k=1;k<=j;k++) s2();
for(int k=1;k<2*n;k++) for(int l=1;l<2*n;l++) myplus(g[k][l],f[k][l]);
}
}
}
int main()
{
freopen("tree.in","r",stdin);
freopen("tree.out","w",stdout);
cin>>n>>m>>q;
for(int i=1;i<=m;i++)
{
int u,v,w;cin>>u>>v>>w;
e[u][v]=e[v][u]=w;
}
pre();det::pre();F::pre();G::pre();
int C=(n-2)*(n-1)/2;
while(q--)
{
LL t;
scanf("%lld",&t);
t%=mod-1;
int ans=0;
for(int i=1;i<2*n;i++)for(int j=1;j<2*n;j++)if(g[i][j])
{
int v=((i-n+mod+C)%mod+1ll*(j-n+mod)*w%mod)%mod;
myplus(ans,1ll*g[i][j]*Pow(v,t)%mod);
}
ans=1ll*ans*Inv(Pow(n*(n-1)/2,t))%mod;
printf("%d\n",ans);
}
return 0;
}
/*
3 2 3
1 2 1
2 3 2
0
1
2
*/