自然数幂和的伯努利数做法是对于指数型母函数的精彩应用。
普通的母函数是因为对于等比数列的研究可以迁移到幂级数上,从而拓展生成函数的运算而简化运算,一般化数列的运算。
指数型母函数是将泰勒展开的研究迁移到幂/阶乘级数(???)而规定生成函数的运算,从而达到描述数列运算的目的。
这个证明很简洁。
B
(
z
)
=
z
e
z
−
1
B(z) = \frac z{e^z-1}
B(z)=ez−1z
也就是
B
(
z
)
e
z
−
1
z
=
1
B(z)\frac {e^z - 1}z = 1
B(z)zez−1=1
∑
i
=
0
n
B
i
i
!
×
1
(
n
−
i
+
1
)
!
=
[
n
=
1
]
\sum_{i=0}^n \frac {B_i}{i!}\times \frac{1}{(n-i+1)!} = [n=1]
∑i=0ni!Bi×(n−i+1)!1=[n=1]
两边同时
×
n
!
\times n!
×n!就是
∑
i
=
0
n
B
i
(
n
+
1
i
)
=
[
n
=
1
]
\sum_{i=0}^n B_i \binom{n+1}i = [n=1]
∑i=0nBi(in+1)=[n=1]
然后 答案就是
B
(
z
)
B(z)
B(z) 与
e
n
z
−
1
z
\frac {e^{nz} - 1}z
zenz−1的卷积的第
k
+
1
k+1
k+1项。
记住
B
(
z
)
=
∑
B
i
x
i
i
!
B(z) = \frac {\sum B_ix ^ i}{i!}
B(z)=i!∑Bixi
e
n
z
−
1
z
=
∑
n
i
x
i
−
1
i
!
\frac {e^{nz} -1}z = \frac {\sum n^ix ^ {i-1}}{i!}
zenz−1=i!∑nixi−1
可以推出
S
k
k
!
=
∑
i
=
0
k
B
i
i
!
∗
n
k
+
1
−
i
(
k
+
1
−
i
)
!
\frac {S_k} {k!} = \sum_{i=0}^k \frac {Bi}{i!} * \frac {n^{k+1-i}} {(k+1-i)!}
k!Sk=∑i=0ki!Bi∗(k+1−i)!nk+1−i
即
S
k
=
1
k
+
1
∑
i
=
0
k
(
k
+
1
i
)
B
i
n
k
+
1
−
i
S_k = \frac 1{k+1}\sum_{i=0}^k \binom{k+1}i B_i n^{k+1-i}
Sk=k+11∑i=0k(ik+1)Bink+1−i
这就是网上的神奇公式。
这里要注意一下,伯努利数不一定是整数,所以这个式子不可以拿来证明1~n的自然数幂和被n整除。
这个式子的好处在于,你只需要计算一次伯努利数,之后都可以用,但是对于模数的要求。。。。
请用第二类斯特林数。
为什么伯努利数分子上有一个
z
z
z?
因为不这样的话就有一项
z
−
1
z^{-1}
z−1
对于伯努利数更自然的描述式应该是这个:
B
n
=
∑
i
=
0
n
B
i
(
n
i
)
,
n
>
1
B_n = \sum_{i=0}^n B_i\binom{n}{i} , n>1
Bn=∑i=0nBi(in),n>1
可以等价于指数生成函数的:
B
(
x
)
+
x
=
B
(
x
)
e
x
B(x) + x = B(x) e^x
B(x)+x=B(x)ex
左边由于
n
>
1
n>1
n>1的限制需要
+
x
+x
+x
那么可以得到
B
(
x
)
=
x
e
x
−
1
B(x) = \frac x{e^x-1}
B(x)=ex−1x
但是递推的时候因为
B
n
B_n
Bn会被消掉,用这个式子算
B
n
−
1
B_{n-1}
Bn−1。。。。
记得
B
0
=
1
B_0 =1
B0=1
51nod 1228 O(k^2)伯努利数模板。
#include<bits/stdc++.h>
#define mod 1000000007
#define maxn 2005
#define LL long long
using namespace std;
int B[maxn],C[maxn][maxn],inv[maxn];
int main()
{
inv[0] = inv[1] = 1;
for(int i=0;i<maxn;i++)
{
if(i > 1) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
C[i][0] = 1;
for(int j=1;j<=i;j++)
C[i][j] = (C[i-1][j-1] + C[i-1][j]) % mod;
}
B[0] = 1;
for(int i=1;i<maxn-1;i++)
{
for(int j=0;j<i;j++) B[i] = (B[i] - 1ll * B[j] * C[i+1][j]) % mod;
B[i] = 1ll * B[i] * inv[i+1] % mod;
}
int T;
for(scanf("%d",&T);T--;)
{
LL n , k;
scanf("%lld%lld",&n,&k);
n %= mod;
int ans = 0;
for(int i=1,pw=n+1;i<=k+1;i++,pw=(1ll * pw * n + pw) % mod)
ans = (ans + 1ll * C[k+1][i] * B[k+1-i] % mod * pw) % mod;
printf("%lld\n",((1ll * ans * inv[k+1]) % mod + mod) % mod);
}
}
51 nod 1258 V4
多项式求逆
n
l
o
g
n
nlogn
nlogn
#include<bits/stdc++.h>
#define maxn 200005
#define M ((1<<15)-1)
#define mod 1000000007
#define LL long long
using namespace std;
const double Pi = 3.1415926535897932384626433832795;
struct cplx
{
double r,i;
cplx(double r=0,double i=0):r(r),i(i){}
cplx operator +(const cplx &B)const{ return cplx(r+B.r,i+B.i); }
cplx operator -(const cplx &B)const{ return cplx(r-B.r,i-B.i); }
cplx operator *(const cplx &B)const{ return cplx(r*B.r-i*B.i,r*B.i+i*B.r); }
cplx conj(){ return cplx(r,-i); }
}w[maxn];
int r[maxn],wlen;
inline void FFT(cplx A[maxn],int lgn,int tp)
{
int n = 1<<lgn;
for(int i=1;i<n;i++) r[i] = (r[i>>1]>>1)|((i&1)<<(lgn-1));
for(int i=1;i<n;i++) if(i<r[i]) swap(A[i],A[r[i]]);
for(int L=2;L<=n;L<<=1)
for(int st=0,l=L>>1,inc=wlen/l;st<n;st+=L)
for(int k=0,loc=0;k<l;k++,loc+=inc)
{
cplx tmp = (tp==1?w[loc]:w[loc].conj())*A[st+k+l];
A[st+k+l]=A[st+k]-tmp,A[st+k]=A[st+k]+tmp;
}
if(tp==-1) for(int i=0;i<n;i++) A[i].r/=n,A[i].i/=n;
}
cplx s[4][maxn];
inline void mul(int a[maxn],int b[maxn],int lgn,int c[maxn])
{
if(lgn == 14)
lgn = 14;
int n = 1<<lgn;
for(int i=0;i<n;i++) s[0][i]=cplx(a[i]>>15,b[i]>>15),s[1][i]=cplx(a[i]&M,b[i]&M);
FFT(s[0],lgn,1),FFT(s[1],lgn,1);
for(int i=0;i<n;i++)
{ cplx a[4]={s[0][i],s[0][(n-i)&(n-1)].conj(),s[1][i],s[1][(n-i)&(n-1)].conj()};
cplx b[4]={(a[0]+a[1])*cplx(0.5,0),(a[0]-a[1])*cplx(0,-0.5),(a[2]+a[3])*cplx(0.5,0),(a[2]-a[3])*cplx(0,-0.5)};
s[2][i]=b[0]*b[1]+cplx(0,1)*b[2]*b[3],s[3][i]=b[0]*b[3]+cplx(0,1)*b[1]*b[2]; }
FFT(s[2],lgn,-1),FFT(s[3],lgn,-1);
for(int i=0;i<n;i++)
{ LL a[4]={llround(s[2][i].r)%mod,llround(s[2][i].i)%mod,llround(s[3][i].r)%mod,llround(s[3][i].i)%mod};
c[i]=((a[0]<<30)+a[1]+((a[2]+a[3])<<15))%mod; }
}
int Pow(int base,int k)
{ int ret=1;for(;k;k>>=1,base=1ll*base*base%mod) if(k&1)ret=1ll*ret*base%mod; return ret; }
void Inv(int a[maxn],int lgn,int b[maxn])
{ if(lgn==0){ b[0]=Pow(a[0],mod-2);return; }
Inv(a,lgn-1,b);
int n = 1<<(lgn+1);
static int tmp[3][maxn];
for(int i=0;i<n;i++)
{ if(i<(n>>1)) tmp[0][i] = a[i]; else tmp[0][i] = 0;
if(i<(n>>2)) tmp[1][i] = b[i]; else tmp[1][i] = 0; }
mul(tmp[0],tmp[1],lgn+1,tmp[2]);
for(int i=0;i<n;i++) tmp[2][i] = (i==0) * 2 - tmp[2][i];
mul(tmp[2],tmp[1],lgn+1,b);
for(int i=(n>>1);i<n;i++) b[i]=0;
}
int a[maxn],B[maxn],inv[maxn]={1,1},invf[maxn]={1,1};
int main()
{
int lgn = 16; wlen = 1<<16;
for(int i=0;i<65536;i++)
{
if(i>1)inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod,
invf[i] = 1ll * invf[i-1] * inv[i] % mod;
w[i] = cplx(cos(i*Pi/wlen),sin(i*Pi/wlen));
}
for(int i=0;i<65536;i++) a[i] = invf[i+1];
Inv(a,lgn,B);
for(int i=0,fac=1;i<65536;i++,fac=1ll*fac*i%mod)
B[i] = (1ll * B[i] * fac % mod + mod)%mod;
int T;
for(scanf("%d",&T);T--;)
{
LL n , k;
scanf("%lld%lld",&n,&k);
n %= mod;
int ans = 0;
for(int i=1,pw=(n+1),C=k+1;i<=k+1;i++,pw=(1ll*pw*n+pw)%mod,C=1ll*C*inv[i]%mod*(k+2-i)%mod)
ans = (ans + 1ll * C * B[k+1-i] % mod * pw) % mod;
printf("%d\n",(1ll*ans*inv[k+1]%mod+mod)%mod);
}
}