ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

WD与积木

2021-04-01 21:01:08  阅读:239  来源: 互联网

标签:WD limits 积木 int dfrac sum hat


XVI.WD与积木

本题有两种思路。

首先,两种思路共同的地方在于都将期望化成了\(\dfrac{\text{所有方案一共的层数}}{\text{总共的方案数}}\)。我们设其为\(\dfrac{f_n}{g_n}\)。

思路1:从DP开始

我们先考虑求出\(g_n\)。

我们有

\[g_n=\sum\limits_{i=1}^n\dbinom{n}{i}g_{n-i} \]

它的意义是枚举第一层里放了哪些积木。

拆开组合数,得到

\[\dfrac{g_n}{n!}=\sum\limits_{i=1}^n\dfrac{1}{i!}\times\dfrac{g_{n-i}}{(n-1)!} \]

我们考虑设\(G_i=\dfrac{g_n}{n!},H_i=\dfrac{1}{i!}\),

则即有

\[G_n=\sum\limits_{i=1}^nH_iG_{n-i} \]

到这里已经可以分治FFT求出了。但是我们还有更好的方法。

我们计算\(GH\),就得到了

\[(GH)_n=\sum\limits_{i=0}^nH_iG_{n-i} \]

带回上面的递推式,我们发现这实际上就是

\[(GH)_n=\sum\limits_{i=1}^nH_iG_{n-i}+H_0G_n=2G_n \]

但是这个递推式有一个特例,即\(n=0\)。此时有\((GH)_0=2G_0-1\)。

故最终有

\[GH=2G-1 \]

处理就得到了

\[G=\dfrac{1}{2-H} \]

求逆即可求出\(G\)。

我们再考虑求出\(f_n\)。

我们有

\[f_n=\sum\limits_{i=1}^n\dbinom{n}{i}(f_{n-i}+g_{n-i}) \]

其中,\(i\)还是枚举第一层里放了多少积木;\(f_{n-i}\)是剔除第一层后剩下部分的总共积木数;\(g_{n-i}\)因为对于剩下\(n-i\)块积木,它们的每一种方法都增加了\(1\)层。

还是仿照我们之前的样子,拆开组合数并设\(F_i=\dfrac{f_i}{i!}\),就有

\[F_n=\sum\limits_{i=1}^nH_i(F_{n-i}+G_{n-i}) \]

计算\(H(F+G)\),有

\[\Big(H(F+G)\Big)_n=\sum\limits_{i=0}^nH_i(F_{n-i}+G_{n-i})=2F_n+G_n \]

照例看一下\(f_0\),发现此时仍有\(\Big(H(F+G)\Big)_0=2F_0+G_0\)(\(F_0=0,G_0=H_0=1\))

故最终有

\[H(F+G)=2F+G \]

于是

\[F=\dfrac{G-GH}{H-2} \]

代入\(G=\dfrac{1}{2-H}\),最终就得到

\[F=G(G-1) \]

然后最终答案就是\(\dfrac{f_n}{g_n}=\dfrac{F_n}{G_n}\)。

此种方案的代码:

#include<bits/stdc++.h>
using namespace std;
const int M=1e5;
const int N=1<<20;
const int mod=998244353;
const int all=18;
int ksm(int x,int y){
	int rt=1;
	for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)rt=(1ll*rt*x)%mod;
	return rt;
}
namespace Poly{
	int rev[N];
	const int G=3;
	void NTT(int *a,int tp,int LG){
		int lim=(1<<LG),invlim=ksm(lim,mod-2);
		for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
		for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int md=1;md<lim;md<<=1){
			int rt=ksm(G,(mod-1)/(md<<1));
			if(tp==-1)rt=ksm(rt,mod-2);
			for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
				int w=1;
				for(int i=0;i<md;i++,w=(1ll*w*rt)%mod){
					int x=a[pos+i],y=(1ll*w*a[pos+md+i])%mod;
					a[pos+i]=(x+y)%mod;
					a[pos+md+i]=(x-y+mod)%mod;
				}
			}
		}
		if(tp==-1)for(int i=0;i<lim;i++)a[i]=(1ll*a[i]*invlim)%mod;
	}
	int A[N],B[N],C[N];
	void mul(int *a,int *b,int *c,int LG){//using: Array A and B
		int lim=(1<<LG);
		for(int i=0;i<lim;i++)A[i]=B[i]=0;
		for(int i=0;i<(lim>>1);i++)A[i]=a[i],B[i]=b[i];
		NTT(A,1,LG),NTT(B,1,LG);
		for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
		NTT(A,-1,LG);
		for(int i=0;i<lim;i++)c[i]=A[i];
	}
	void inv(int *a,int *b,int LG){//using: Array C
		b[0]=ksm(a[0],mod-2);
		for(int k=1;k<=LG+1;k++){
			mul(b,a,C,k);
			for(int i=0;i<(1<<k);i++)C[i]=(mod-C[i])%mod;
			(C[0]+=2)%=mod;
			mul(C,b,b,k);
		}
	}	
}
int a[N];
int f[N],g[N],h[N],res[N];
int fac[N],inv[N];
int T,n;
int main(){
	fac[0]=1;
	for(int i=1;i<=M;i++)fac[i]=1ll*fac[i-1]*i%mod;
	inv[M]=ksm(fac[M],mod-2);
	for(int i=M-1;i>=0;i--)inv[i]=1ll*inv[i+1]*(i+1)%mod;
	for(int i=0;i<=M;i++)h[i]=inv[i];
	for(int i=0;i<=M;i++)a[i]=mod-h[i];
	(a[0]+=2)%=mod;
	Poly::inv(a,g,all);
	for(int i=0;i<=M;i++)a[i]=g[i];
	(a[0]+=mod-1)%=mod;
	Poly::mul(a,g,f,all+1);
	for(int i=0;i<=M;i++)res[i]=1ll*f[i]*ksm(g[i],mod-2)%mod;
	scanf("%d",&T);
	while(T--)scanf("%d",&n),printf("%d\n",res[n]);
	return 0;
}

思路2:从生成函数开始。

仍然先考虑\(g_n\)。

我们考虑一个一层放多少积木的EGF,\(F(x)=\sum\limits_{i=1}^n\dfrac{1}{i!}=e^x-1\)

EGF的卷积就是将这些积木标号的过程;所以我们最终的答案就是

\[\hat{g}=\sum\limits_{i=0}^{\infty}F^i \]

它是\(F\)的一阶前缀和,则有

\[\hat{g}=\dfrac{1}{1-F}=\dfrac{1}{2-e^x} \]

通过复原\(e^x\)为\(\sum\limits_{i=0}^n\dfrac{1}{i!}\)可以直接通过多项式求逆求出\(\hat{g}\)。

我们再考虑\(f_n\)。

则我们只需要乘上一个层数即可;则有

\[\hat{f}=\sum\limits_{i=0}^{\infty}iF^i \]

它是\(F\)的二阶前缀和;于是就有

\[\hat{f}=\dfrac{F}{(1-F)^2}=\dfrac{e^x-1}{(2-e^x)^2} \]

直接求出即可。

所以我们最终要求的就是\(\dfrac{f_n}{g_n}=\dfrac{\hat{f}_n}{\hat{g}_n}\)。

代码就不贴了。

标签:WD,limits,积木,int,dfrac,sum,hat
来源: https://www.cnblogs.com/Troverld/p/14608054.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有