【题意】给定n，求Σi=0~nΣj=1~i s(i,j)*2^j*j!，n<=10^5。
【算法】生成函数+排列组合+多项式求逆
【题解】参考： [BZOJ4555][Tjoi2016&Heoi2016]求和-NTT-多项式求逆
$ans=\sum_{i=0}^{n}\sum_{j=0}^{i}s(i,j)*2^j*j!$
令$g(n)=\sum_{j=0}^{n}s(n,j)*2^j*j!$
则ans是Σg(i)，只要计算出g(i)的生成函数就可以统计答案。
g(n)可以理解为将n个数划分成若干集合，每个集合有2个属性的排列数。基于此实际意义，通过枚举第一个集合来递推g(n)。
$g(n)=\sum_{i=1}^{n}2*C(n,i)*g(n-i)$
特别的，g(0)=1。
两边乘n!（令人窒息的操作），得
$\frac{g(n)}{n!}=\sum_{i=1}^{n}\frac{2}{i!}*\frac{g(n-i)}{(n-i)!}$
这已经是卷积的形式了：
$F(n)=\sum_{i=1}^{n}\frac{2}{n!}*x^i$
$G(n)=\sum_{i=0}^{n}\frac{g(n)}{n!}*x^i$
注意此时F*G卷积后，G(0)的位置是0，所以
$G(n)=F(n)G(n)+1$
移项得，
$G(n)=\frac{1}{1-F(n)}$
最后，多项式求逆即可，复杂度O(n log n)。
#include<cstdio> #include<cstring> #include<algorithm> using namespace std; const int maxn=300010,MOD=998244353; void gcd(int a,int b,int &x,int &y){ if(!b){x=1;y=0;}else{gcd(b,a%b,y,x);y-=x*(a/b);} } int inv(int a){int x,y;gcd(a,MOD,x,y);return (x%MOD+MOD)%MOD;} int power(int x,int k){ int ans=1; while(k){ if(k&1)ans=1ll*ans*x%MOD; x=1ll*x*x%MOD; k>>=1; } return ans; } namespace ntt{ int o[maxn],oi[maxn],f[maxn]; void init(int n){ int x=1,p=power(3,(MOD-1)/n); for(int k=0;k<n;k++){ o[k]=x;oi[k]=inv(o[k]); x=1ll*x*p%MOD; } } void transform(int *a,int n,int *o){ int k=0; while((1<<k)<n)k++; for(int i=0;i<n;i++){ int t=0; for(int j=0;j<k;j++)if((1<<j)&i)t|=(1<<(k-j-1)); if(i<t)swap(a[i],a[t]); } for(int l=2;l<=n;l*=2){ int m=l>>1; for(int *p=a;p!=a+n;p+=l){ for(int i=0;i<m;i++){ int t=1ll*p[i+m]*o[n/l*i]%MOD; p[i+m]=(p[i]-t+MOD)%MOD; p[i]=(p[i]+t)%MOD; } } } } void dft(int *a,int n){transform(a,n,o);} void idft(int *a,int n){ transform(a,n,oi);int v=inv(n); for(int i=0;i<n;i++)a[i]=1ll*a[i]*v%MOD; } void pinv(int *F,int *g,int n){ if(n==1){g[0]=inv(F[0]);return;} pinv(F,g,n>>1);n<<=1; init(n);// for(int i=0;i<n/2;i++)f[i]=F[i],f[i+n/2]=0; dft(f,n);dft(g,n); for(int i=0;i<n;i++)g[i]=1ll*g[i]*(2-1ll*f[i]*g[i]%MOD+MOD)%MOD;//1ll idft(g,n);for(int i=n/2;i<n;i++)g[i]=0;// } } int F[maxn],G[maxn],n,fac[maxn]; int main(){ int n,N=1; scanf("%d",&n);n++; while(N<n)N*=2; fac[0]=1; for(int i=1;i<n;i++)fac[i]=1ll*fac[i-1]*i%MOD; for(int i=1;i<n;i++)F[i]=((-2*inv(fac[i]))%MOD+MOD)%MOD; F[0]++; ntt::pinv(F,G,N); int ans=0; for(int i=0;i<n;i++)ans=(ans+1ll*G[i]*fac[i]%MOD)%MOD; printf("%d",ans); return 0; }
注意：多项式求逆过程中每次都要对不同的n进行一次预处理omega[]。
另一种做法
【算法】斯特林数+NTT
【题解】首先有第二类斯特林数的通项公式。
$s(n,m)=\frac{1}{m!}\sum_{k=0}^{m}(-1)^k*C(m,k)*(m-k)^n$
当斯特林数s(n,m)满足m>n时，上述公式计算结果为0，所以第二个Σ可以扩展到n。
$ans=\sum_{i=0}^{n}\sum_{j=0}^{n}s(i,j)*s^j*j!$
代入第二类斯特林数公式。
$ans=\sum_{i=0}^{n}\sum_{j=0}^{n}2^j*j!*\frac{1}{j!}\sum_{k=0}^{j}(-1)^k*\frac{j!}{k!(j-k)!}*(j-k)^i$
通过组合数的分解，向卷积靠拢。
注意到Σi只对最后一个括号有影响，所以移动到最后。
$ans=\sum_{j=0}^{n}2^j*j!\sum_{k=0}^{j}\frac{(-1)^k}{k!}*\frac{\sum_{i=0}^{n}(j-k)^i}{(j-k)!}$
这已经是标准的卷积形式（第二个函数分子是等比数列可以直接计算）。
使用NTT计算。
注意：
1.n=0，只有使0^0=1，斯特林数通项公式才能处理s(0,0)的情况。
2.n=1，等比数列求和公式不能处理Σ1^i即q=1的情况。
所以特殊地，g[0]=1，g[1]=n+1。先计算完再n++。
#include<cstdio> #include<cstring> #include<algorithm> using namespace std; const int maxn=300010,MOD=998244353; int power(int x,int k){ int ans=1; while(k){ if(k&1)ans=1ll*ans*x%MOD; x=1ll*x*x%MOD; k>>=1; } return ans; } int inv(int x){return power(x,MOD-2);} namespace ntt{ int o[maxn],oi[maxn]; void init(int n){ int x=1,g=power(3,(MOD-1)/n); for(int k=0;k<n;k++){ o[k]=x;oi[k]=inv(o[k]); x=1ll*x*g%MOD; } } void transform(int *a,int n,int *o){ int k=0; while((1<<k)<n)k++; for(int i=0;i<n;i++){ int t=0; for(int j=0;j<k;j++)if((1<<j)&i)t|=(1<<(k-j-1)); if(i<t)swap(a[i],a[t]); } for(int l=2;l<=n;l*=2){ int m=l>>1; for(int *p=a;p!=a+n;p+=l){ for(int i=0;i<m;i++){ int t=1ll*p[i+m]*o[n/l*i]%MOD; p[i+m]=(p[i]-t+MOD)%MOD; p[i]=(p[i]+t)%MOD; } } } } void dft(int *a,int n){transform(a,n,o);} void idft(int *a,int n){ transform(a,n,oi); int x=inv(n); for(int i=0;i<n;i++)a[i]=1ll*a[i]*x%MOD; } } int n,fac[maxn],f[maxn],g[maxn]; int main(){ scanf("%d",&n); fac[0]=1; for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%MOD; for(int i=0;i<=n;i++)f[i]=1ll*((i&1)?MOD-1:1)*inv(fac[i])%MOD; for(int i=2;i<=n;i++)g[i]=1ll*(1-power(i,n+1)+MOD)*inv((1-i+MOD)%MOD)%MOD*inv(fac[i])%MOD; g[0]=1;g[1]=n+1;// n++;int N=1;// while(N<n+n)N*=2; ntt::init(N); ntt::dft(f,N);ntt::dft(g,N); for(int i=0;i<N;i++)f[i]=1ll*f[i]*g[i]%MOD; ntt::idft(f,N); int ans=0; for(int i=0;i<n;i++)ans=(ans+1ll*power(2,i)*fac[i]%MOD*f[i]%MOD)%MOD; printf("%d",ans); return 0; }
