题解P3711：【仓鼠的数学题】

这题黑的丫！怎么会掉紫呢！

noteskey

伯努利数... 这里 有介绍哟~ 写的非常详细呢~

反正这题就是推柿子...

另外就是黈力算法的运用 QWQ

我们令 \(ANS(x)\) 为答案多项式，那么这个多项式可以这么求：

（下面我们定义 \(S(n,k)\) 为自然幂和函数（不是第二类斯特林数！），即 \(\sum_{i=0}^{ n} i^k\)）

\[\begin{aligned}ANS(x)=& \sum_{k=0}^{n} S_k (x) a_k \\ =& \sum_{k=0}^{n} (S(x,k)+x^k) a_k \\ =& \sum_{k=0}^n a_k \Big( x^k +{1\over k+1} \sum_{i=0}^k \begin{pmatrix} k+1\\i \end{pmatrix} B_i x^{k+1-i}\Big) \\ =& \sum_{k=0}^n a_k \Big( x^k + k! \sum_{i=0}^k {B_i\over i!} {x^{k+1-i} \over (k+1-i)! }\Big) \\ =& \sum_{k=0}^n a_k x^k + \sum_{k=0}^n a_k k! \sum_{i=0}^k {B_i\over i!} {x^{k+1-i} \over (k+1-i)! } \\ =& \sum_{i=0}^n a_i x^i + \sum_{i=1}^{n+1} x^i i! \sum_{k=i-1}^{n} a_k k!{B_{k+1-i}\over (k+1-i)!} \end{aligned} \]

然后感觉做不下去了呢...后面虽说像是卷积的形式然鹅根本就不是卷积呢QWQ

首先，我们看着表达式太长了，于是乎把表达式换成一个函数：

\[\begin{aligned}f(i)=&~a_ii! \\g(i) =&~{B_i \over i!} \end{aligned} \]

那么原来的式子就是：

\[ANS=\sum_{i=0}^n a_i x^i + \sum_{i=1}^{n+1} x^i i! \sum_{k=i-1}^{n} f(k)g(k+1-i) \]

这样的话，我们就更加清晰的发现后面的式子不是卷积的形式了 XD

那么我们令 \(gr(i)=g(n-i)\) （也就是翻转多项式）

然后原来的式子就是：

\[ANS=\sum_{i=0}^n a_i x^i+ \sum_{i=1}^{n+1} x^i i! \sum_{k=i-1}^{n} f(k)gr(n+i-k-1) \]

发现这里就是最高项改成 \(n+i\) 项的卷积了，于是我们发现 要计算后面的 \(\sum_{k=i-1}^n\) 什么的算出来的是 \(f* gr\) 的第 \((n+i-1)\) 项

于是我们把 \(f* gr\) 的多项式求出来，拿第 \((n-i+1)\) 项乘上 \(i!\) ，再加上 \(a_i\) 就是最后多项式的第 i 项答案了

好像没有黑题难度？ 前提你得知道伯努利数这玩意儿丫！ 【雾

而且这个推导过程...好像比 shadowice 巨巨的短很多丫...【小声

code

不知道为什么不开 Ofast 跑得比开了快... 600ms +

//by Judge
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define Rg register
#define fp(i,a,b) for(Rg int i=(a),I=(b)+1;i<I;++i)
#define ll long long
using namespace std;
const int mod=998244353;
const int iG=332748118;
const int M=1<<19|3;
typedef int arr[M];
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
char buf[1<<21],*p1=buf,*p2=buf;
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int dec(int x,int y){return x<y?x-y+mod:x-y;}
inline int inc(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int read(){ int x=0,f=1; char c=getchar();
	for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
	for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} char sr[1<<21],z[20];int CCF=-1,Z;
inline void Ot(){fwrite(sr,1,CCF+1,stdout),CCF=-1;}
inline void print(int x,char chr=' '){
    if(CCF>1<<20)Ot();if(x<0)sr[++CCF]=45,x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++CCF]=z[Z],--Z);sr[++CCF]=chr;
} int n,limit; arr fac,ifac,a,B,C,r;
inline int qpow(int x,int p=mod-2,int s=1){
	for(;p;p>>=1,x=mul(x,x)) if(p&1) s=mul(s,x); return s;
}
inline void init(int n){ int len=-1;
	for(limit=1;limit<=n;limit<<=1) ++len;
	fp(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<len); 
}
inline void NTT(int* a,int tp){ fp(i,0,limit-1) if(i<r[i]) swap(a[i],a[r[i]]);
	for(Rg int mid=1;mid<limit;mid<<=1){ Rg int I=mid<<1,Gn=qpow(tp?3:iG,(mod-1)/I);
		for(Rg int j=0,x;j<limit;j+=I) for(Rg int k=0,g=1;k<mid;++k,g=mul(g,Gn))
			x=mul(a[j+k+mid],g),a[j+k+mid]=dec(a[j+k],x),a[j+k]=inc(a[j+k],x);
	} if(tp) return ; int inv=qpow(limit); fp(i,0,limit-1) a[i]=mul(a[i],inv);
}
void get_inv(int* a,int* b,int n){
	if(n==1) return b[0]=qpow(a[0]),void(); static int c[M],d[M];
	get_inv(a,b,n>>1),init(n); for(int i=0;i<n;++i) c[i]=a[i],d[i]=b[i];
	for(int i=n;i<limit;++i) c[i]=d[i]=0; NTT(c,1),NTT(d,1);
	fp(i,0,limit-1) c[i]=mul(c[i],mul(d[i],d[i])); NTT(c,0);
	fp(i,0,n-1) b[i]=dec(inc(b[i],b[i]),c[i]);
}
inline void prep(int len){ B[0]=ifac[0]=ifac[1]=fac[0]=fac[1]=1;
	fp(i,2,len) ifac[i]=mul(mod-mod/i,ifac[mod%i]);
	fp(i,2,len) ifac[i]=mul(ifac[i-1],ifac[i]),fac[i]=mul(fac[i-1],i);
	get_inv(ifac+1,B,len);
}
int main(){ n=read();
	for(limit=1;limit<=n;limit<<=1); prep(limit);
	fp(i,0,n) a[i]=read(),C[i]=mul(a[i],fac[i]);
	reverse(B,B+1+n),init(n<<1|1);
	fp(i,n+1,limit-1) B[i]=C[i]=0; NTT(B,1),NTT(C,1);
	fp(i,0,limit-1) B[i]=mul(B[i],C[i]); NTT(B,0);
	fp(i,0,n+1) B[n+i-1]=mul(B[n+i-1],ifac[i]);
	fp(i,0,n) B[n+i-1]=inc(B[n+i-1],a[i]);
	print(a[0]); fp(i,1,n+1) print(B[n+i-1]);
	return sr[CCF]='\n',Ot(),0;
}

more...

其实我们只需要用另外一种伯努利数，就可以让推导更加简洁（没简洁多少，代码复杂度也基本没变的说）

这个伯努利数就是之前提到的那篇博客里面讲的 \(B^+\)，这个伯努利数列满足的性质更加适合做这道题...

于是我们令 \(B=B^+\)， 然后重推一遍：

（其实没必要，给个代码就 OJBK 了？ 【滑稽）

\[\begin{aligned}ANS(x)=& \sum_{k=0}^{n} S_k (x) a_k \\ =& \sum_{k=0}^{n} S^+(x,k) a_k \\ =& \sum_{k=0}^n a_k \Big( {1\over k+1} \sum_{i=0}^k \begin{pmatrix} k+1\\i \end{pmatrix} B_i x^{k+1-i} \Big) \\ =& \sum_{k=0}^n a_k k! \sum_{i=0}^k {B_i\over i!} {x^{k+1-i} \over (k+1-i)! } \\ =& \sum_{i=1}^{n+1} x^i i! \sum_{k=i-1}^{n} a_k k!{B_{k+1-i}\over (k+1-i)!} \end{aligned} \]

依然令：

\[\begin{aligned}f(i)=&~a_ii! \\g(i) =&~{B_i \over i!} \\g_r(i)=&~g(n-i)\end{aligned} \]

（注意这里的 \(B_i\) 就是 \(B_i^+\) ）

\[ANS=\sum_{i=1}^{n+1} x^i i! \sum_{k=i-1}^{n} f(k)gr(n+i-k-1) \]

同样的，把 \(f* g_r\) 的第 \((n-i+1)\) 项乘上 \(i!\) 就是最后的答案了 ，不同的就是这里不需要再加上 \(a_i\) XD

code*2

代码...可谓没有什么区别

//by Judge
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define Rg register
#define fp(i,a,b) for(Rg int i=(a),I=(b)+1;i<I;++i)
#define ll long long
using namespace std;
const int mod=998244353;
const int iG=332748118;
const int M=1<<19|3;
typedef int arr[M];
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
char buf[1<<21],*p1=buf,*p2=buf;
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int dec(int x,int y){return x<y?x-y+mod:x-y;}
inline int inc(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int read(){ int x=0,f=1; char c=getchar();
	for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
	for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} char sr[1<<21],z[20];int CCF=-1,Z;
inline void Ot(){fwrite(sr,1,CCF+1,stdout),CCF=-1;}
inline void print(int x,char chr=' '){
    if(CCF>1<<20)Ot();if(x<0)sr[++CCF]=45,x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++CCF]=z[Z],--Z);sr[++CCF]=chr;
} int n,limit; arr fac,ifac,a,B,C,r;
inline int qpow(int x,int p=mod-2,int s=1){
	for(;p;p>>=1,x=mul(x,x)) if(p&1) s=mul(s,x); return s;
}
inline void init(int n){ int len=-1;
	for(limit=1;limit<=n;limit<<=1) ++len;
	fp(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<len); 
}
inline void NTT(int* a,int tp){ fp(i,0,limit-1) if(i<r[i]) swap(a[i],a[r[i]]);
	for(Rg int mid=1;mid<limit;mid<<=1){ Rg int I=mid<<1,Gn=qpow(tp?3:iG,(mod-1)/I);
		for(Rg int j=0,x;j<limit;j+=I) for(Rg int k=0,g=1;k<mid;++k,g=mul(g,Gn))
			x=mul(a[j+k+mid],g),a[j+k+mid]=dec(a[j+k],x),a[j+k]=inc(a[j+k],x);
	} if(tp) return ; int inv=qpow(limit); fp(i,0,limit-1) a[i]=mul(a[i],inv);
}
void get_inv(int* a,int* b,int n){
	if(n==1) return b[0]=qpow(a[0]),void(); static int c[M],d[M];
	get_inv(a,b,n>>1),init(n); for(int i=0;i<n;++i) c[i]=a[i],d[i]=b[i];
	for(int i=n;i<limit;++i) c[i]=d[i]=0; NTT(c,1),NTT(d,1);
	fp(i,0,limit-1) c[i]=mul(c[i],mul(d[i],d[i])); NTT(c,0);
	fp(i,0,n-1) b[i]=dec(inc(b[i],b[i]),c[i]);
}
inline void prep(int len){
	B[0]=ifac[0]=ifac[1]=fac[0]=fac[1]=1;
	fp(i,2,len) ifac[i]=mul(mod-mod/i,ifac[mod%i]);
	fp(i,2,len) ifac[i]=mul(ifac[i-1],ifac[i]);
	fp(i,2,len) fac[i]=mul(fac[i-1],i);
	get_inv(ifac+1,B,len); B[1]=mod-B[1];
}
int main(){ n=read();
	for(limit=1;limit<=n;limit<<=1); prep(limit);
	fp(i,0,n) a[i]=read(),C[i]=mul(a[i],fac[i]);
	reverse(B,B+1+n),init(n<<1|1);
	fp(i,n+1,limit-1) B[i]=C[i]=0; NTT(B,1),NTT(C,1);
	fp(i,0,limit-1) B[i]=mul(B[i],C[i]); NTT(B,0);
	fp(i,0,n+1) B[n+i-1]=mul(B[n+i-1],ifac[i]);
	print(a[0]); fp(i,1,n+1) print(B[n+i-1]);
	return sr[CCF]='\n',Ot(),0;
}