[P5162] WD与积木
每种堆法(理解成名次序列,举例3,3,8,2和7,7,100,2都对应2,2,1,3这个名次序列)等概率出现;题目中“两种堆法不同当且仅当某个积木在两种堆法中处于不同的层中”可见这是个组合问题,于是设一个名次的生成函数f(x)=sum i=1->inf 1*x^i/(i!) 。因为f(x)的常数项为0,有f^inf(x)=0。
名次序列的方案数的生成函数:a(x) = sum i=0->inf f^i(x) = 1/(1-f(x))
名次序列的各种方案的名次个数之和的生成函数:b(x) = sum i=0->inf i*f^i(x) = f(x)/(1-f(x))^2 = a^2(x)-a(x) = a(x)*(a(x)-1)
预处理出a(x),b(x)的前max_n项系数,对于n,答案为b(x)[n]/a(x)[n]。
注意事项
关于从0开始求和
如果是从1开始的情形:
a(x)=sum i=1->inf f^i(x)
= f(x)*(1-f^inf(x))/(1-f(x))
= f(x)/(1-f(x))
= (1-(1-f(x)))/(1-f(x))
= 1/(1-f(x))-1b(x)=sum i=1-> inf i*f^i(x)
= ...
= f(x)/(1-f(x))^2-1抛开a(x)[0]和b(x)[0],似乎关系并没有什么影响。
关于egf的卷积运算
某个生成函数f(x)=sum i=0->inf a[i]*x^i/(i!)
对于egf意义下的卷积
f(x)*f(x)=sum i=0->inf (sum j=0->i c(i,j)*a[j]*a[i-j]) x^i/(i!)
= sum i=0->inf (sum j=0->i (i!)/(j!)/((i-j)!)*a[j]*a[i-j]) x^i/(i!)
= sum i=0->inf (sum j=0->i a[j]/(j!)*a[i-j]/((i-j)!)) x^i
与在ogf意义下的卷积
f(x)=sum i=0->inf b[i]*x^i ,b[i]=a[i]/(i!)
f(x)*f(x)=sum i=0->inf (sum j=0->i b[j]*b[i-j]) x^i
完全一致.这告诉我们在使用egf的具体计算时,可以将1/(i!)放到系数中去用一般多项式算法计算得到一个ogf结构,
所得到的函数的各系数再乘上i!就是正确的egf结构下的标志序列了;显然,该性质也适用于egf的连续卷积!换句话说,egf的标志序列尽管除以一个i!,然后与ogf无异地各种操作,最终结果再让标志序列乘上i!就好了。
代码实现中因为除法关系没有乘上i!。
参考实现
#include <bits/stdc++.h> using namespace std; const int max_n=1e5+10; const int p=998244353,g=3; int a[max_n<<2]; int b[max_n<<2]; int c[max_n<<2]; int f[max_n<<2]; int r[max_n<<2]; int frr[max_n]; int qpow(long long x,int y) { long long c=1; for(; y; y>>=1, x=x*x%p) if(y&1) c=c*x%p; return c; } void ntt(int *a,int lmt,int tp) { for(int i=0; i<lmt; ++i) if(i<r[i]) swap(a[i],a[r[i]]); for(int m=1; m<lmt; m<<=1) { int wm=qpow(g,(p-1)/(m<<1)); if(tp==-1) wm=qpow(wm,p-2); for(int i=0; i<lmt; i+=(m<<1)) { long long w=1, tmp; for(int j=0; j<m; ++j, w=w*wm%p) { tmp=w*a[i+j+m]%p; a[i+j+m]=(a[i+j]-tmp+p)%p; a[i+j]=(a[i+j]+tmp)%p; } } } if(tp==-1) { long long tmp=qpow(lmt,p-2); for(int i=0; i<lmt; ++i) a[i]=tmp*a[i]%p; } } void polyrev(int deg,int *a,int *b) { if(deg==1) { b[0]=qpow(a[0],p-2); return; } polyrev((deg+1)>>1,a,b); int lmt=1, l=0; while(lmt<(deg<<1)) lmt<<=1, ++l; for(int i=0; i<lmt; ++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0; i<deg; ++i) c[i]=a[i]; fill(c+deg,c+lmt,0); ntt(c,lmt,1), ntt(b,lmt,1); for(int i=0; i<lmt; ++i) b[i]=(2ll-1ll*c[i]*b[i]%p+p)%p*b[i]%p; ntt(b,lmt,-1); fill(b+deg,b+lmt,0); } int main() { long long fac=1; for(int i=1; i<max_n; ++i) fac=fac*i%p; frr[max_n-1]=qpow(fac,p-2); for(int i=max_n-1; i; --i) frr[i-1]=1ll*i*frr[i]%p; for(int i=0; i<max_n; ++i) a[i]=p-frr[i]; a[0]=(a[0]+2)%p; polyrev(max_n,a,b); for(int i=0; i<max_n; ++i) a[i]=f[i]=b[i]; int lmt=1; while(lmt<=max_n) lmt<<=1; lmt<<=1; a[0]=(a[0]+p-1)%p; ntt(f,lmt,1), ntt(a,lmt,1); for(int i=0; i<lmt; ++i) f[i]=1ll*f[i]*a[i]%p; ntt(f,lmt,-1); int t,n; scanf("%d",&t); while(t--) { scanf("%d",&n); printf("%lld\n",1ll*f[n]*qpow(b[n],p-2)%p); } return 0; }
感谢@weng_weiji的答疑,思路也是来自于他(她?)的题解;代码参考了@owen_codeisking。
也感谢来自其它大佬的恶意嘲讽
上一篇: 台风用伞
下一篇: 中国最大餐饮机器人制造商将进军日本