1.定义式
定义伯努利数列(B_n)满足：
2.递推式
可以发现定义式里面包含了(B_n)这一项，于是把(B_n)提出来：
直接用定义式求是(O(n^2))的复杂度
3.生成函数
把定义式的循环上界减一，得：
注意到组合数上标变成了(n)，再加个(B_n)：
组合数很烦，把它拆开来：
两边都写成生成函数的形式：
设伯努利数的指数型生成函数为(B(x)=sum_{n=0}frac{B_n}{n!}x^n)，那么左边显然就是(B(x)e^x)，右边就是(B(x))。
但是细想却不对劲，如果(B(x)e^x=B(x))，那么(e^x=1)，显然不成立。注意到一开始把循环上标减了1，而根据定义，伯努利数定义式的(n)必须大于(0)，所以上面的式子只有当(n-1>0)时成立。
于是考虑列出(B(x)e^x)的第(0,1)项：
而(B(x))的第一项只有(B_1x)，少了个(B_0x)，加上去即可：
可以解出(B(x)=frac{x}{e^x-1})。
注意到分母的第(0)次项为(0)，所以把分子除下去之后，第(0)次项不为(0)，可以用多项式求逆(O(nlogn))做。
4.应用
说了这么多，伯努利数有什么用？
答：可以求解自然数幂和。
定义自然数幂和(S(m,k)=sum_{i=1}^mi^k)，直接求单个是(O(nlogn))的（快速幂是(O(logn))）。使用伯努利数可以对于固定的(m)做到(O(nlogn))求出前(n)个。
下面说下求法：
对于固定的(m)，还是设(S(m,k))的指数型生成函数为
交换求和符号得到：
注意到分母和(B(x))一样，把(B(x))带进去：
把(e^{mx})展开得到：
前面这个(B(x)e^x)不好处理，但注意到(B(x)e^x=B(x)+x)，于是引入新的(B'(x)=B(x)+x)，把(B'(x))展开：
注意到(S_m(x)=sum_{k=0}S(m,k)frac{x^k}{k!})，所以：
这个就可以称为自然数幂和的通项公式。
注意到通项公式可以写成卷积的形式：
因为(B'(x)=B(x)+x)，可以(O(nlogn))求出每一项(B'_n)，所以可以(O(nlogn))求出(m)固定时(kin[1,n])的(S(m,k))值
5.题目
luoguP3711仓鼠的数学题
已知(a_0...a_n)，设(S_k(x)=sum_{i=0}^xi^k)，求出(sum_{k=0}^nS_k(x)a_k)，并输出其每一项的系数，可以证明答案是个(n+1)次多项式。(nleq 250000)，(0^0=1)，答案对(998244353)取模。
可以发现这里的(S_k(x))相比上面所说的多了个下标(0)。但是发现当(k>0)时，(0^k=0)，所以(S_k(x)=sum_{i=1}^xi^k(k>0))；而当(k=0)时，(0^k=1)，所以(S_k(x)=sum_{i=1}^xi^k+1(k=0))。
所以可以重新定义(S_k(x))为(sum_{i=1}^xi^k)，这样答案的式子就是：
只需要关注前面的求和式，用伯努利数将其展开得到：
那么设(f_i=a_ii!)，(g_i=frac{B'_i}{i!})，引入(g'_i=g_{n-i})，那么上面的式子可以接着写成：
注意到(g_i=frac{B'_i}{i!}=frac{i![x^i]B'(x)}{i!}=[x^i]B'(x)=[x^i](frac{x}{e^x-1}+x))，可以(O(nlogn))求，答案的卷积也是(O(nlogn))，卷积的第(n+i-1)次项除以(i!)就是答案的第(i)项。
注意到上面式子没有第(0)次项，所以答案第(0)次项就是(a_0)。
#include<bits/stdc++.h>
#define rg register
#define il inline
#define cn const
#define gc getchar()
#define fp(i,a,b) for(rg int i=(a),ed=(b);i<=ed;++i)
#define fb(i,a,b) for(rg int i=(a),ed=(b);i>=ed;--i)
#define add(a,b) (((a)+(b))%mod)
#define inc(a,b) (((a)-(b)+mod)%mod)
#define mul(a,b) (1ll*(a)*(b)%mod)
#define div(a,b) (1ll*(a)*fpow((b,mod-2))%mod)
using namespace std;
typedef cn int cint;
il int rd(){
rg int x(0); rg char c(gc);
while(c<'0'||'9'<c)c=gc;
while('0'<=c&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=gc;
return x;
}
cint maxn=2000010,mod=998244353,G=3,invG=332748118;
int n,lim=1,a[maxn],fac[maxn],ifac[maxn],bo[maxn],f[maxn],g[maxn];
il int fpow(int a,int b,int ans=1){
for(;b;b>>=1,a=mul(a,a)){
if(b&1)ans=mul(ans,a);
}
return ans;
}
namespace poly{
int A[maxn],B[maxn];
int lim,l,rev,r[maxn];
il void init(cint &l1,cint &l2){
lim=1,l=0;
while(lim<=l1+l2)lim<<=1,++l;
rev=fpow(lim,mod-2);
fp(i,0,lim-1)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
il void ntt(int *a,cint &f){
fp(i,0,lim-1)if(i<r[i])swap(a[i],a[r[i]]);
for(rg int md=1;md<lim;md<<=1){
rg int len=md<<1,Gn=fpow(f?G:invG,(mod-1)/len);
for(rg int l=0;l<lim;l+=len){
rg int Pow=1;
for(rg int nw=0;nw<md;++nw,Pow=mul(Pow,Gn)){
rg int x=a[l+nw],y=mul(Pow,a[l+nw+md]);
a[l+nw]=add(x,y),a[l+nw+md]=inc(x,y);
}
}
}
}
void getinv(int *a,int *f,int n){
if(n==1){
f[0]=fpow(a[0],mod-2);
return;
}
getinv(a,f,n>>1);
init(n-1,n-1);
fp(i,0,n-1)A[i]=a[i],B[i]=f[i];
fp(i,n,lim)A[i]=B[i]=0;
ntt(A,1),ntt(B,1);
fp(i,0,lim)A[i]=mul(B[i],inc(2,mul(A[i],B[i])));
ntt(A,0);
fp(i,0,n-1)f[i]=mul(A[i],rev);
}
il void mult(int *a,int *b,int *c,cint &l1,cint &l2,cint &l3){
init(l1,l2);
fp(i,0,l1)A[i]=a[i];
fp(i,0,l2)B[i]=b[i];
fp(i,l1+1,lim)A[i]=0;
fp(i,l2+1,lim)B[i]=0;
ntt(A,1),ntt(B,1);
fp(i,0,lim)A[i]=mul(A[i],B[i]);
ntt(A,0);
fp(i,0,l3)c[i]=mul(A[i],rev);
}
}
int main(){
n=rd();
fp(i,0,n)a[i]=rd();
fac[0]=1;
fp(i,1,n+1)fac[i]=mul(fac[i-1],i);
ifac[n+1]=fpow(fac[n+1],mod-2);
fb(i,n+1,1)ifac[i-1]=mul(ifac[i],i);
fp(i,0,n)g[i]=ifac[i+1];
while(lim<=n)lim<<=1;
poly::getinv(g,bo,lim),++bo[1];
fp(i,0,n)f[i]=mul(a[i],fac[i]);
fp(i,0,n)g[i]=bo[n-i];
poly::mult(f,g,f,n,n,n<<1);
printf("%d",a[0]);
fp(i,1,n+1)printf(" %lld",mul(f[n+i-1],ifac[i]));
return 0;
}