luogu P3711 仓鼠的数学题

背景：

发现自己推公式的能力好差。

题目传送门：

https://www.luogu.org/problemnew/show/P3711

题意：

求${\sum }_{k=0}^n{a}_k{\sum }_{i=0}^x{i}^k$的每一项的系数。
你不需要管$x$是什么，毒瘤出题人就是让你不不知道某些知识读不懂题。

思路：

大力推式子。
$$\sum _{k=0}^n{a}_k\sum _{i=0}^x{i}^k$$
考虑伯努利数，带入${\sum }_{i=0}^x{i}^k=0^k+x^k+{\sum }_{i=1}^{x-1}{i}^k=x^k+{\sum }_{i=1}^{x-1}=x^k+\frac{1}{k+1}{\sum }_{i=1}^{k+1}{C}_{k+1}^i\cdot B_{k+1-i}\cdot (x-1+1{)}^i=x^k+\frac{1}{k+1}{\sum }_{i=1}^{k+1}{C}_{k+1}^i\cdot B_{k+1-i}\cdot x^i$有：
$$=\sum _{k=0}^n{a}_k(x^k+\frac{1}{k+1}\sum _{i=1}^{k+1}{C}_{k+1}^i\cdot B_{k+1-i}\cdot x^i)$$

考虑$C$是有对称性的，因此有：
$$\begin{array}{rl}& =\sum _{k=0}^n{a}_k(x^k+\frac{1}{k+1}\sum _{i=0}^k{C}_{k+1}^i\cdot B_i\cdot x^{k+1-i})\\ & =\sum _{k=0}^n{a}_k(x^k+\frac{1}{k+1}\sum _{i=0}^k\frac{(k+1)!}{i!(k+1-i)!}\cdot B_i\cdot x^{k+1-i})\end{array}$$

抵消$k+1$，得：
$$\begin{array}{rl}& =\sum _{k=0}^n{a}_k(x^k+\sum _{i=1}^{k+1}\frac{k!}{i!(k+1-i)!}\cdot B_i\cdot x^{k+1-i})\\ & =\sum _{k=0}^n{a}_k{x}^k+\sum _{k=0}^n\sum _{i=0}^k{a}_k\frac{k!}{i!(k+1-i)!}\cdot B_i\cdot x^{k+1-i}\end{array}$$

设$d=k+1-i$，有：
$$=\sum _{d=0}^n{a}_d{x}^d+\sum _{d=1}^{n+1}\sum _{k=d-1}^n{a}_k\frac{k!}{(k+1-d)!d!}\cdot B_{k+1-d}\cdot x^d$$

设$f_i=a_i\cdot i!,g=\frac{B_{k+1-d}}{(k+1-d)!}$所以有：
$$\begin{array}{rl}& =\sum _{d=1}^{n+1}{a}_d{x}^d+\sum _{d=0}^n\sum _{k=d-1}^n\frac{f_k{g}_{k+1-d}}{d!}\cdot x^d\\ & =\sum _{d=1}^{n+1}{a}_d{x}^d+\sum _{d=0}^n\frac{1}{d!}\cdot x^d\sum _{k=d-1}^n{f}_k{g}_{k+1-d}\end{array}$$

设$h_i=g_{n-i}$，则有：$g_{k+1-d}=h_{n-(k+1-d)}=h_{n-k+d-1}$，即：
$$=\sum _{d=1}^{n+1}{a}_d{x}^d+\sum _{d=0}^n\frac{1}{d!}\cdot x^d\sum _{k=d-1}^n{f}_k{h}_{n-k+d-1}$$

发现后面是一个卷积的形式，最高项的次数为$k+(n-k+d-1)=n+d-1$，用第$n+d-1$的系数乘上$\frac{1}{d!}$，最后加上$a_k$即可。
特别地：$an{s}_0=a_0$（$0^0=1$）。

代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
typedef long double LD;
const int mod=998244353,G=3,inv_G=332748118;
	int r[2000010],A[2000010],a[2000010],b[2000010],f[2000010],g[2000010];
	int fac[2000010],Inv[2000010],B[2000010];
	int limit,n,l;
	LL sum;
int dg(int x,int k)
{
	if(!k) return 1;
	int op=dg(x,(k>>1));
	if(k&1) return (LL)op*op%mod*x%mod; else return (LL)op*op%mod;
}
int inv(int x)
{
	return dg(x,mod-2);
}
void NTT(int *now,int limit,int op)
{
	for(int i=0;i<limit;i++)
		if(i<r[i]) swap(now[i],now[r[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		int wn=dg(op==1?G:inv_G,(mod-1)/(mid<<1));
		for(int j=0;j<limit;j+=(mid<<1))
		{
			int w=1;
			for(int k=0;k<mid;k++,w=(int)(((LL)w*wn)%mod))
			{
				LL x=now[j+k],y=(LL)w*now[j+k+mid]%mod;
				now[j+k]=((LL)x+y)%mod;
				now[j+k+mid]=((LL)x-y+mod)%mod;
			}
		}
	}
}
void init(int n)
{
	limit=1,l=0;
	while(limit<(n<<1))
		limit<<=1,l++;
	for(int i=1;i<limit;i++)
		r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
}
void dft(int *f,int n,int limit)
{
	NTT(f,limit,-1);
	int INV=inv(limit);
	for(int i=0;i<limit;i++)
		f[i]=(LL)f[i]*INV%mod;
}
void poly_inv(int *f,int *g,int n)
{
	if(n==1)
	{
		g[0]=inv(f[0]);
		return;
	}
	poly_inv(f,g,(n+1)>>1);
	init(n);
	memset(a,0,sizeof(a));
	memset(b,0,sizeof(b));
	for(int i=0;i<n;i++)
		a[i]=f[i],b[i]=g[i];
	NTT(a,limit,1),NTT(b,limit,1);
	for(int i=0;i<limit;i++)
		b[i]=(LL)b[i]*((2ll-(LL)a[i]*b[i]%mod+mod)%mod)%mod;
	dft(b,n,limit);
	for(int i=0;i<n;i++)
		g[i]=b[i];
}
void INIT(int n)
{
	fac[0]=fac[1]=1;
	Inv[0]=Inv[1]=1;
	for(int i=2;i<=n;i++)
	{
		fac[i]=(LL)fac[i-1]*i%mod;
		Inv[i]=((LL)mod-mod/i)*Inv[mod%i]%mod;
	}
	for(int i=2;i<=n;i++)
		Inv[i]=(LL)Inv[i-1]*Inv[i]%mod;
	for(int i=0;i<=n;i++)
		f[i]=Inv[i+1];
	poly_inv(f,g,n);
	for(int i=0;i<=n;i++)
		B[i]=(LL)g[i]*fac[i]%mod;
}
int calc_C(int n,int m)
{
	return (LL)fac[n]*Inv[n-m]%mod*Inv[m]%mod;
}
int main()
{
	scanf("%d",&n);
	for(int i=0;i<=n;i++)
		scanf("%d",&A[i]);
	init(n);
	INIT(limit);
	for(int i=0;i<=n;i++)
	{
		f[i]=(LL)A[i]*fac[i]%mod;
		g[n-i]=(LL)B[i]*Inv[i]%mod;
	}
	init(n+1);
	memset(a,0,sizeof(a));
	memset(b,0,sizeof(b));
	for(int i=0;i<=n;i++)
		a[i]=f[i],b[i]=g[i];
	NTT(a,limit,1),NTT(b,limit,1);
	for(int i=0;i<limit;i++)
		a[i]=(LL)a[i]*b[i]%mod;
	dft(a,n+1,limit);
	for(int i=0;i<=limit;i++)
		g[i]=a[i];
	printf("%d",A[0]);
	for(int d=1;d<=n+1;d++)
		printf(" %d",((LL)A[d]+(LL)Inv[d]*g[n+d-1]%mod)%mod);
}