伯努利数与自然数幂和

自然数幂和的伯努利数做法是对于指数型母函数的精彩应用。
普通的母函数是因为对于等比数列的研究可以迁移到幂级数上,从而拓展生成函数的运算而简化运算,一般化数列的运算。
指数型母函数是将泰勒展开的研究迁移到幂/阶乘级数(???)而规定生成函数的运算,从而达到描述数列运算的目的。
在这里插入图片描述
这个证明很简洁。
B ( z ) = z e z − 1 B(z) = \frac z{e^z-1} B(z)=ez1z
也就是 B ( z ) e z − 1 z = 1 B(z)\frac {e^z - 1}z = 1 B(z)zez1=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×(ni+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 zenz1的卷积的第 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!} zenz1=i!nixi1
可以推出
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+1i)!nk+1i
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+11i=0k(ik+1)Bink+1i
这就是网上的神奇公式。
这里要注意一下,伯努利数不一定是整数,所以这个式子不可以拿来证明1~n的自然数幂和被n整除。
这个式子的好处在于,你只需要计算一次伯努利数,之后都可以用,但是对于模数的要求。。。。
请用第二类斯特林数。

为什么伯努利数分子上有一个 z z z?
因为不这样的话就有一项 z − 1 z^{-1} z1
对于伯努利数更自然的描述式应该是这个:
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)=ex1x

但是递推的时候因为 B n B_n Bn会被消掉,用这个式子算 B n − 1 B_{n-1} Bn1。。。。
记得 B 0 = 1 B_0 =1 B0=1

关于 B + B^+ B+ B − B^- B

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);
	}
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值