伯努利数学习笔记 & 51nod 1228 & 51nod 1258

置顶 zsyz_ZZY

于 2019-07-06 15:07:57 发布

阅读量291

点赞数 1

分类专栏： # NTT # 倍增 # 伯努利数文章标签：数论 NTT 倍增

本文链接：https://blog.csdn.net/zsyz_ZZY/article/details/94858270

版权

NTT 同时被 3 个专栏收录

20 篇文章 2 订阅

订阅专栏

倍增

18 篇文章 0 订阅

订阅专栏

伯努利数

2 篇文章 1 订阅

订阅专栏

背景：

最近做到了一些 $\sum_{i=1}^{n}i^k$ 的题目，学习了伯努利数这种操作。

题目传送门：

https://www.51nod.com/Challenge/Problem.html#!#problemId=1228
https://www.51nod.com/Challenge/Problem.html#!#problemId=1258

题意：

多组询问求： $\sum_{i=1}^{n}i^k$ 。

正题：

定义 $B$ 表示伯努利数。
$B_0=1,B_1=-\frac{1}{2}，B_2=\frac{1}{6},B_3=0,B_4=-\frac{1}{30}$

$B_5=0,B_6=\frac{1}{42}，B_7=0,B_8=-\frac{1}{30},B_9=0$

$\text{Part}1$ ：

其实伯努利数有 $n^2$ 的求法。
$B_n=[n=0]-\frac{1}{n+1}\sum_{i=0}^{n-1}B_iC_{n+1}^{i}$

~~证明不会啊，好像也没有什么用吧。~~
递推即可。

$\text{Part}2$ ：

伯努利数的生成函数的定义是：
$\frac{x}{e^x-1}=\sum_{n=0}^{\infty}B_n\frac{x^n}{n!}$
其中 $x \in C$ ， $B_i$ 为伯努利数的第 $i$ 项。
因为： $\sum_{i=0}^{\infty}\frac{x^i}{i!}=e^x$ ，所以有：
$\frac{x}{e^x-1}=\frac{x}{-1+\sum_{i=0}^{\infty}\frac{x^i}{i!}}$

当 $i = 0$ 时， $\frac{x_i}{i!}=\frac{1}{1}=1$ ，所以有：
$=\frac{x}{\sum_{i=1}^{\infty}\frac{x^i}{i!}}$

上面的 $x$ 相当于整体右移一位，得：
$=\frac{1}{\sum_{i=0}^{\infty}\frac{x^i}{(i+1)!}}$

因此可以用 $\text{NTT}$ 多项式求逆。

$\text{Part}3$ ：

伯努利数能解决什么问题呢？
有一条这样的公式。
$\sum_{i=1}^{n}i^{k}=\frac{1}{k+1}\sum_{i=1}^{k+1}C_{k+1}^{i}\cdot B_{k+1-i}\cdot (n+1)^i$
另一条长这样：
$\sum_{i=1}^{n}i^{k}=\frac{1}{k+1}\sum_{i=0}^{k}C_{k+1}^{i}\cdot B_{i}\cdot (n+1)^{k+1-i}$
假设我们求出了伯努利数，我们就可以 $\Theta(k)$ 求解了。
~~时间复杂度极其优秀（然而预处理需要更大的时间复杂度）。~~

代码 $1$ ：

$\text{P1228}$ ：
$\Theta(n^2)$ 的：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
#define mod 1000000007
using namespace std;
	LL n;
	int k;
	LL ans;
	LL fac[2010],inv[2010],Inv[2010],B[2010];
LL calc_C(LL n,LL m)
{
	if(m>n) return 0;
	return n<mod?fac[n]*Inv[n-m]%mod*Inv[m]%mod:calc_C(n/mod,m/mod)*calc_C(n%mod,m%mod)%mod;
}
void init(LL ma)
{
	fac[0]=fac[1]=inv[0]=inv[1]=Inv[0]=Inv[1]=1;
	for(int i=2;i<=ma;i++)
		fac[i]=fac[i-1]*i%mod;
	for(int i=2;i<=ma+1;i++) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
	for(int i=2;i<=ma;i++) Inv[i]=Inv[i-1]*inv[i]%mod;

	B[0]=1;
	for(int i=1;i<=ma;i++)
	{
		for(int j=0;j<i;j++)
			B[i]=(B[i]+B[j]*calc_C(i+1,j)%mod)%mod;
		B[i]=(-inv[i+1]*B[i]%mod+mod)%mod;
	}
}
int main()
{
	int T;
	init(2005);
	scanf("%d",&T);
	while(T--)
	{
		scanf("%lld %d",&n,&k);
		n%=mod;
		ans=0;
		LL op=1;
		for(int i=1;i<=k+1;i++)
		{
			op=op*(n+1)%mod;
			ans=(ans+calc_C(k+1,i)*B[k+1-i]%mod*op%mod)%mod;
		}
		printf("%lld\n",ans*inv[k+1]%mod);
	}	
}

代码 $2$ ：

$\text{P1258}$ ：
$\Theta(n\log n)$ 的：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define R register
#define I inline
using namespace std;
typedef long long LL;
typedef long double LD;
const int Mod[3]={469762049,998244353,1004535809},MOD=1000000007,G=3;
LL M=(LL)Mod[0]*Mod[1];
	int a[300010],b[300010],f[300010],g[300010],ans[3][300010];
	int fac[300010],Inv[300010],B[300010];
	int limit,n,l,r[300010];
	int mod,inv_G;
	LL sum;
I int read()
{
	int x=0,f=1;
	char ch=getchar();
	for(;ch<'0'||ch>'9';ch=getchar())
		if(ch=='-') f=-1;
	for(;ch>='0'&&ch<='9';x=(x<<3)+(x<<1)+(ch^48),ch=getchar());
	return x*f;
}
LL dg(int x,int k,LL mod)
{
	if(!k) return 1;
	LL op=dg(x,(k>>1),mod);
	if(k&1) return op*op%mod*x%mod; else return op*op%mod;
}
I LL times(LL a,LL b,LL mod)
{
	a%=mod,b%=mod;
	return ((a*b-(LL)((LL)((LD)a/mod*b+1e-3)*mod))%mod+mod)%mod;
}
I int inv(int x,int mod)
{
	return dg(x,mod-2,(LL)mod);
}
I void NTT(int *now,int limit,int op)
{
	for(R int i=0;i<limit;i++)
		if(i<r[i]) swap(now[i],now[r[i]]);
	for(R int mid=1;mid<limit;mid<<=1)
	{
		int wn=(int)dg(op==1?G:inv_G,(mod-1)/(mid<<1),mod);
		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=(int)((LL)w*now[j+k+mid]%mod);
				now[j+k]=(int)(((LL)x+y)%mod);
				now[j+k+mid]=(int)(((LL)x-y+mod)%mod);
			}
		}
	}
}
I int merge(LL x,LL y,LL z)
{
	LL A=(times(x*Mod[1]%M,dg(Mod[1]%Mod[0],Mod[0]-2,Mod[0]),M)+
		  times(y*Mod[0]%M,dg(Mod[0]%Mod[1],Mod[1]-2,Mod[1]),M))%M;
	LL K=((z-A)%Mod[2]+Mod[2])%Mod[2]*dg(M%Mod[2],Mod[2]-2,Mod[2])%Mod[2];
	return (int)(((K%MOD)*(M%MOD)%MOD+A%MOD)%MOD);
}
I void mul(int *f,int *g,int *h,int n,int limit,int md)
{
	for(int i=n;i<limit;i++)
		a[i]=b[i]=h[i]=0;
	for(R int i=0;i<n;i++)
		a[i]=f[i],b[i]=g[i];
	mod=md;
	inv_G=inv(G,mod);
	NTT(a,limit,1);
	NTT(b,limit,1);
	for(R int i=0;i<limit;i++)
		a[i]=((LL)a[i]*b[i])%mod;
	NTT(a,limit,-1);
	int INV=(int)inv(limit,mod);
	for(R int i=0;i<n;i++)
		h[i]=(LL)a[i]*INV%mod;
}
void work(int n,int *f,int *g)
{
	if(n==1)
	{
		g[0]=inv(f[0],MOD);
		return;
	}
	work((n+1)>>1,f,g);
	limit=1,l=0;
	while(limit<(n<<1))
		limit<<=1,l++;
	for(R int i=1;i<limit;i++)
		r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
	mul(f,g,ans[0],n,limit,Mod[0]),mul(f,g,ans[1],n,limit,Mod[1]),mul(f,g,ans[2],n,limit,Mod[2]);
	ans[0][0]=ans[1][0]=ans[2][0]=1;
	for(R int i=1;i<=n;i++)
		ans[0][i]=ans[1][i]=ans[2][i]=(-merge(ans[0][i],ans[1][i],ans[2][i])+MOD)%MOD;
	mul(ans[0],g,ans[0],n,limit,Mod[0]),mul(ans[1],g,ans[1],n,limit,Mod[1]),mul(ans[2],g,ans[2],n,limit,Mod[2]);
	for(R int i=0;i<n;i++)
		g[i]=merge(ans[0][i],ans[1][i],ans[2][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];
	work(n,f,g);
	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()
{
	int T,y;
	LL x;
	limit=1;
	init(50000+1);
	scanf("%d",&T);
	while(T--)
	{
		sum=0;
		scanf("%lld %d",&x,&y);
		x%=MOD;
		int op=1;
		for(int i=1;i<=y+1;i++)
		{
			op=(LL)op*(x+1)%MOD;
			sum=((LL)sum+(LL)calc_C(y+1,i)*B[y+1-i]%MOD*op%MOD)%MOD;
		}
		printf("%d\n",(LL)sum*inv(y+1,MOD)%MOD);
	}
}

zsyz_ZZY

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
伯努利数学习笔记 & 51nod 1228 & 51nod 1258

背景：最近做到了一些∑i=1nik\sum_{i=1}^{n}i^k∑i=1nik的题目，学习了伯努利数这种操作。题目传送门：https://www.51nod.com/Challenge/Problem.html#!#problemId=1228题意：多组询问求：∑i=1nik\sum_{i=1}^{n}i^k∑i=1nik。正题：定义BBB表示伯努利数。B0=1...
复制链接

扫一扫