首页 > 代码库 > hdu5824 graph
hdu5824 graph
传送门
题意：定义一个无向图的权值为图中形为树的连通块数量的$k$次方，求所有$n$个点有标号的简单无向图的权值之和。
这个题还是很妙的啊……(好吧，其实只有最后的复合函数求导比较有意思……)
先套路一发，定义答案的$EGF$为$F(x)$一棵树的$EGF$为$T(x)$，每个连通块都不是树的$EGF$为$G(x)$，强制连通的无向图的$EGF$为$H(x)$，显然有
\begin{align}F(x)=\left(\sum_{i\ge 0}\frac{i^k T(x)^i}{i!}\right)G(x)\end{align}
意思就是枚举树的数量，统计数量之后乘上贡献，左边要除以$i!$是为了去重(因为几个树组成的图是集合，不是序列，元素没有顺序)。
根据一些简单的知识不难得到以下结论：
$[x^n]T(x)=n^{n-2}$
$G(x)=\exp(H(x)-T(x))$(因为每个连通块都不能是树，那么我们用连通图个数减掉树的个数，再用这个东西搞一个集合即可，不难发现就是对一个连通块的$EGF$做一下$\exp$的结果。)
$H(x)$可以跑多项式$\ln$之类的东西得到(参见城市规划的做法)，那么后面的$G(x)$就好算了，然后再解决前面的那个东西就可以$O(n)$得出最后的答案了(因为只要求一项，所以暴力求出卷积的第$n$项即可)。
定义
\begin{align}A_k(x)=\sum_{i\ge 0}\frac{i^k T(x)^i}{i!}\end{align}
注意到$A_k(x)$其实是一个复合函数，那么假设有$f(T)=A_k(x)$，我们可以尝试对$f(T)$求个导，而根据复合函数求导法则($(f(g(x)))’=f’(g(x))g’(x)$)有$[T^i]f’(T)=\frac{i^k T(x)^{i-1}T‘(x)i}{i!}=T‘(x)\frac{i^{k+1}T(x)^{i-1}}{i!}$，因此$[T^i]f’(T)=\frac{i^k T(x)^{i-1}T‘(x)i}{i!}=T‘(x)\frac{i^{k+1}T(x)^{i-1}}{i!}$。
因为$f(T)$本质就是$A_k(x)$，因此$f’(T)$和$A_k’(x)$是等价的，所以有
\begin{align}A_k‘(x)=\frac{T‘(x)}{T(x)}\sum_{i\ge 1}\frac{i^{k+1}T(x)^i}{i!}=\frac{T‘(x)}{T(x)}A_{k+1}(x)\end{align}
(注意这里不用处理常数项的问题，因为$k>0$时有$[x^0]A_k(x)=0$)
也就是说$A_{k+1}(x)=A_k‘(x)\frac{T(x)}{T‘(x)}$，那么直接利用这个递推$k$次即可，边界也不难得到：
\begin{align}A_0(x)=\sum_{i\ge 0}\frac{i^0 T(x)^i}{i!}=\exp(T(x))\end{align}
然后就可以做了，复杂度$O(kn\log n)$，因为$\exp$只需要做两次，所以常数还是不大的。
(虽然比起NTT优化DP的做法来说代码长还跑得慢……)
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 using namespace std; 5 const int maxn=65550,p=998244353,g=3; 6 void NTT(int*,int,int); 7 void getexp(int*,int*,int); 8 void getln(int*,int*,int); 9 void getinv(int*,int*,int); 10 void getderivative(int*,int*,int); 11 void getintegrate(int*,int*,int); 12 int qpow(int,int,int); 13 int n,N=1,k,A[maxn],B[maxn],C[maxn],T[maxn],G[maxn],fac[maxn],fac_inv[maxn],inv[maxn],ans=0; 14 int main(){ 15 scanf("%d%d",&n,&k); 16 while(N<=n)N<<=1; 17 fac[0]=fac_inv[0]=1; 18 for(int i=1;i<N;i++)fac[i]=(long long)fac[i-1]*i%p; 19 fac_inv[N-1]=qpow(fac[N-1],p-2,p); 20 for(int i=N-2;i;i--)fac_inv[i]=(long long)fac_inv[i+1]*(i+1)%p; 21 for(int i=1;i<N;i++)inv[i]=(long long)fac_inv[i]*fac[i-1]%p; 22 T[1]=G[0]=1; 23 for(int i=2;i<=n;i++)T[i]=(long long)qpow(i,i-2,p)*fac_inv[i]%p; 24 for(int i=1;i<=n;i++)G[i]=(long long)qpow(2,(((long long)i*(i-1))>>1)%(p-1),p)*fac_inv[i]%p; 25 getln(G,B,N); 26 for(int i=0;i<N;i++)B[i]=(B[i]-T[i]+p)%p; 27 getexp(B,G,N); 28 fill(B,B+(N<<1),0); 29 getderivative(T,B,N); 30 getinv(B,C,N); 31 copy(T,T+N,B); 32 NTT(B,N<<1,1); 33 NTT(C,N<<1,1); 34 for(int i=0;i<(N<<1);i++)C[i]=(long long)B[i]*C[i]%p; 35 NTT(C,N<<1,-1); 36 fill(C+N,C+(N<<1),0); 37 NTT(C,N<<1,1); 38 getexp(T,A,N); 39 for(int j=1;j<=k;j++){ 40 fill(B,B+(N<<1),0); 41 getderivative(A,B,N); 42 NTT(B,N<<1,1); 43 for(int i=0;i<(N<<1);i++)B[i]=(long long)B[i]*C[i]%p; 44 NTT(B,N<<1,-1); 45 copy(B,B+N,A); 46 } 47 for(int i=0;i<=n;i++)ans=(ans+(long long)A[i]*G[n-i]%p)%p; 48 printf("%d",(int)((long long)ans*fac[n]%p)); 49 return 0; 50 } 51 void NTT(int *A,int n,int tp){ 52 for(int i=1,j=0,k;i<n-1;i++){ 53 k=n; 54 do j^=(k>>=1);while(j<k); 55 if(i<j)swap(A[i],A[j]); 56 } 57 for(int k=2;k<=n;k<<=1){ 58 int wn=qpow(g,(tp>0?(p-1)/k:(p-1)/k*(long long)(p-2)%(p-1)),p); 59 for(int i=0;i<n;i+=k){ 60 int w=1; 61 for(int j=0;j<(k>>1);j++,w=(long long)w*wn%p){ 62 int a=A[i+j],b=(long long)w*A[i+j+(k>>1)]%p; 63 A[i+j]=(a+b)%p; 64 A[i+j+(k>>1)]=(a-b+p)%p; 65 } 66 } 67 } 68 if(tp<0){ 69 int inv=qpow(n,p-2,p); 70 for(int i=0;i<n;i++)A[i]=(long long)A[i]*inv%p; 71 } 72 } 73 void getexp(int *A,int *C,int n){ 74 static int B[maxn]; 75 fill(C,C+(n<<1),0); 76 C[0]=1; 77 for(int k=2;k<=n;k<<=1){ 78 getln(C,B,k); 79 for(int i=0;i<k;i++)B[i]=(A[i]-B[i]+p)%p; 80 B[0]=(B[0]+1)%p; 81 NTT(B,k<<1,1); 82 NTT(C,k<<1,1); 83 for(int i=0;i<(k<<1);i++)C[i]=(long long)B[i]*C[i]%p; 84 NTT(C,k<<1,-1); 85 fill(C+k,C+(k<<1),0); 86 } 87 } 88 void getln(int *A,int *C,int n){ 89 static int B[maxn]; 90 getderivative(A,B,n); 91 fill(B+n,B+(n<<1),0); 92 getinv(A,C,n); 93 NTT(B,n<<1,1); 94 NTT(C,n<<1,1); 95 for(int i=0;i<(n<<1);i++)B[i]=(long long)B[i]*C[i]%p; 96 NTT(B,n<<1,-1); 97 getintegrate(B,C,n); 98 fill(C+n,C+(n<<1),0); 99 } 100 void getinv(int *A,int *C,int n){ 101 static int B[maxn]; 102 fill(C,C+(n<<1),0); 103 C[0]=qpow(A[0],p-2,p); 104 for(int k=2;k<=n;k<<=1){ 105 copy(A,A+k,B); 106 fill(B+k,B+(k<<1),0); 107 NTT(B,k<<1,1); 108 NTT(C,k<<1,1); 109 for(int i=0;i<(k<<1);i++)C[i]=C[i]*((2-(long long)C[i]*B[i]%p+p)%p)%p; 110 NTT(C,k<<1,-1); 111 fill(C+k,C+(k<<1),0); 112 } 113 } 114 void getderivative(int *A,int *C,int n){ 115 for(int i=1;i<n;i++)C[i-1]=(long long)A[i]*i%p; 116 C[n-1]=0; 117 } 118 void getintegrate(int *A,int *C,int n){ 119 for(int i=1;i<n;i++)C[i]=(long long)A[i-1]*inv[i]%p; 120 C[0]=0; 121 } 122 int qpow(int a,int b,int p){ 123 int ans=1; 124 for(;b;b>>=1,a=(long long)a*a%p)if(b&1)ans=(long long)ans*a%p; 125 return ans; 126 }
ps：代码里把每次乘的$\frac{T(x)}{T’(x)}$算出来之后直接存它的$DFT$了，这样每次乘的时候就只需要对$A_k(x)$做$DFT$和$IDFT$了，可以减小很多常数。
这个题给我的两点启发：
1.看见跟前面那半类似的式子就去试试求导和积分
2.牢记复合函数求导法则(论不学文化课的危害……)
hdu5824 graph