zoukankan      html  css  js  c++  java
  • 伯努利数

    伯努利数

    (B_0=1,B_1=-frac{1}{2},B_2=frac{1}{6},B_3=0,B_4=frac{1}{30})

    可以利用下面的式子计算。

    [B_0=1,sum_{i=0}^nB_iC_{n+1}^i=0 ]

    转化:

    [egin{aligned} sum_{i=0}^nB_iC_{n+1}^i&=0\ sum_{i=0}^{n-1}B_iC_{n}^i&=0\ sum_{i=0}^{n-1}B_iC_{n}^i+B_n&=B_n(n>1)\ sum_{i=0}^{n}B_iC_n^i&=B_n(n>1)\ sum_{i=0}^nfrac{n!}{i!(n-i)!}B_i&=B_n(n>1)\ sum_{i=0}^nfrac{1}{i!(n-i)!}B_i&=frac{B_n}{n!}(n>1)\ end{aligned} ]

    对于(B_i)构建指数型生成函数$$B(x)=sum_{i=0}^{infty}frac{B_i}{i!}x^i$$

    发现上面转化出来的左侧是(B(x))乘上({1})的指数型生成函数。

    所以我们可以得到式子$$B(x)e^x=B(x)+x$$

    需要额外加(x)的原因是上面转化式子中的右侧没有(B_1)

    所以我们可以得到伯努利数的另外一种形式(B(x)=frac{x}{e^x-1})

    所以我们利用多项式求逆在(O(nlogn))计算伯努利数。

    当然,这个式子也可以拆开,然后(O(n^2))计算。

    [egin{aligned} sum_{i=0}^nB_iC_{n+1}^i&=0\ sum_{i=0}^{n-1}B_iC_{n+1}^i+B_nC_{n+1}^n&=0\ -B_n(n+1)&=sum_{i=0}^{n-1}B_iC_{n+1}^i\ B_n&=-frac{1}{n+1}sum_{i=0}^{n-1}B_iC_{n+1}^i end{aligned} ]

    伯努利数可以用来计算自然数的幂和。

    [egin{aligned} S_k(n)&=sum_{i=0}^{n-1}i^k\ &=frac{1}{k+1}sum_{i=0}^kC_{k+1}^iB_in^{k+1-i} end{aligned} ]

    这个式子是这么来的

    [egin{aligned} A(x)&=sum_{i=0}^{infty}(sum_{j=0}^{n-1}j^i)frac{x^i}{i!}=sum_{j=0}^{n-1}sum_{i=0}^{infty}j^ifrac{x^i}{i!}\ &=sum_{j=0}^{n-1}e^{jx}=frac{e^{nx}-1}{e^x-1}\ &=B(x)frac{e^{nx}-1}{x} end{aligned} ]

    然后再拆掉这个式子就可以得到上面的结果(我说其实我并不会拆这个式子)


    以下是(51Nod1228)(O(n^2))预处理伯努利数

    #include<iostream>
    #include<cstdio>
    #include<cstdlib>
    #include<cstring>
    #include<cmath>
    #include<algorithm>
    #include<set>
    #include<map>
    #include<vector>
    #include<queue>
    using namespace std;
    #define ll long long
    #define RG register
    #define MOD 1000000007
    #define MAX 2500
    inline ll read()
    {
        RG ll x=0,t=1;RG char ch=getchar();
        while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
        if(ch=='-')t=-1,ch=getchar();
        while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
        return x*t;
    }
    int fpow(int a,int b)
    {
    	int s=1;
    	while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}
    	return s;
    }
    int B[MAX],jc[MAX],jv[MAX],inv[MAX];
    int C(int n,int m){return 1ll*jc[n]*jv[m]%MOD*jv[n-m]%MOD;}
    int main()
    {
    	B[0]=jc[0]=jv[0]=inv[0]=inv[1]=1;
    	for(int i=2;i<MAX;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
    	for(int i=1;i<MAX;++i)jc[i]=1ll*jc[i-1]*i%MOD;
    	for(int i=1;i<MAX;++i)jv[i]=1ll*jv[i-1]*inv[i]%MOD;
    	for(int i=1;i<MAX;++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;B[i]=(MOD-B[i])%MOD;
    	}
    	int T=read();
    	while(T--)
    	{
    		ll n=read();int k=read(),ans=0;
    		int nw=(n+1)%MOD,q=(n+1)%MOD;
    		for(int i=k;~i;--i,nw=1ll*nw*q%MOD)ans=(ans+1ll*C(k+1,i)*B[i]%MOD*nw)%MOD;
    		ans=1ll*ans*inv[k+1]%MOD;
    		printf("%d
    ",ans);
    	}
    	return 0;
    }
    
    

    一下是(51Nod1258),多项式求逆预处理伯努利数

    #include<iostream>
    #include<cstdio>
    #include<cstdlib>
    #include<cstring>
    #include<cmath>
    #include<algorithm>
    #include<set>
    #include<map>
    #include<vector>
    #include<queue>
    using namespace std;
    #define ll long long
    #define RG register
    #define MOD 1000000007
    #define MAX 150000
    const int NN=50000;
    const int M=sqrt(MOD);
    const double Pi=acos(-1);
    inline ll read()
    {
        RG ll x=0,t=1;RG char ch=getchar();
        while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
        if(ch=='-')t=-1,ch=getchar();
        while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
        return x*t;
    }
    int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
    struct Complex{double a,b;}W[MAX],A1[MAX],A2[MAX],B1[MAX],B2[MAX],X[MAX],Y[MAX],Z[MAX];
    Complex operator+(Complex a,Complex b){return (Complex){a.a+b.a,a.b+b.b};}
    Complex operator-(Complex a,Complex b){return (Complex){a.a-b.a,a.b-b.b};}
    Complex operator*(Complex a,Complex b){return (Complex){a.a*b.a-a.b*b.b,a.b*b.a+a.a*b.b};}
    int r[MAX],N,l;
    void FFT(Complex *P,int N,int opt)
    {
    	for(int i=0;i<N;++i)if(i<r[i])swap(P[i],P[r[i]]);
    	for(int i=1;i<N;i<<=1)
    		for(int p=i<<1,j=0;j<N;j+=p)
    			for(int k=0;k<i;++k)
    			{
    				Complex w=(Complex){W[N/i*k].a,W[N/i*k].b*opt};
    				Complex X=P[j+k],Y=P[i+j+k]*w;
    				P[j+k]=X+Y;P[i+j+k]=X-Y;
    			}
    	if(opt==-1)for(int i=0;i<N;++i)P[i].a/=N;
    }
    void MTT(int *a,int *b,int len,int *c)
    {
    	memset(A1,0,sizeof(A1));memset(B1,0,sizeof(B1));
    	memset(A2,0,sizeof(A2));memset(B2,0,sizeof(B2));
    	for(N=1,l=0;N<=len;N<<=1)++l;
    	for(int i=0;i<N;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    	for(int i=1;i<N;i<<=1)
    		for(int k=0;k<i;++k)W[N/i*k]=(Complex){cos(k*Pi/i),sin(k*Pi/i)};
    	for(int i=0;i<len;++i)a[i]%=MOD,b[i]%=MOD;
    	for(int i=0;i<len;++i)A1[i].a=a[i]/M,A2[i].a=a[i]%M;
    	for(int i=0;i<len;++i)B1[i].a=b[i]/M,B2[i].a=b[i]%M;
    	memset(X,0,sizeof(X));memset(Y,0,sizeof(Y));memset(Z,0,sizeof(Z));
    	FFT(A1,N,1);FFT(A2,N,1);FFT(B1,N,1);FFT(B2,N,1);
    	for(int i=0;i<N;++i)
    	{
    		X[i]=A1[i]*B1[i];
    		Y[i]=A1[i]*B2[i]+A2[i]*B1[i];
    		Z[i]=A2[i]*B2[i];
    	}
    	FFT(X,N,-1);FFT(Y,N,-1);FFT(Z,N,-1);
    	for(int i=0;i<len;++i)
    	{
    		int ans=0;
    		ans=(ll)(X[i].a+0.5)%MOD*M%MOD*M%MOD;
    		ans=(ans+(ll)(Y[i].a+0.5)%MOD*M)%MOD;
    		ans=(ans+(ll)(Z[i].a+0.5))%MOD;
    		c[i]=ans;
    	}
    }
    int c[MAX],d[MAX];
    void Inv(int *a,int *b,int len)
    {
    	if(len==1){b[0]=fpow(a[0],MOD-2);return;}
    	Inv(a,b,len>>1);
    	MTT(a,b,len,c);MTT(b,c,len,d);
    	for(int i=0;i<len;++i)b[i]=(b[i]+b[i])%MOD;
    	for(int i=0;i<len;++i)b[i]=(b[i]+MOD-d[i])%MOD;
    }
    int a[MAX];
    int n,B[MAX],jc[MAX],jv[MAX],inv[MAX];
    int C(int n,int m){return 1ll*jc[n]*jv[m]%MOD*jv[n-m]%MOD;}
    int main()
    {
    	B[0]=jc[0]=jv[0]=inv[0]=inv[1]=1;
    	for(int i=2;i<=NN+NN;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
    	for(int i=1;i<=NN+NN;++i)jc[i]=1ll*jc[i-1]*i%MOD;
    	for(int i=1;i<=NN+NN;++i)jv[i]=1ll*jv[i-1]*inv[i]%MOD;
    	for(int i=0;i<=NN;++i)a[i]=jv[i+1];
    	Inv(a,B,1<<16);
    	for(int i=0;i<=NN;++i)B[i]=1ll*B[i]*jc[i]%MOD;
    	int T=read();
    	while(T--)
    	{
    		int n=read()%MOD,k=read(),ans=0,nw=n+1;
    		for(int i=k;~i;--i,nw=1ll*nw*(n+1)%MOD)ans=(ans+1ll*C(k+1,i)*B[i]%MOD*nw)%MOD;
    		ans=1ll*ans*inv[k+1]%MOD;
    		printf("%d
    ",ans);
    	}
    	return 0;
    }
    
    
  • 相关阅读:
    十大开源Web应用安全测试工具
    HC大会,华为联合合作伙伴发布一站式物联网IoT开发工具小熊派BearPi
    漫谈边缘计算(四):赢家是软还是硬
    漫谈边缘计算(三):5G的好拍档
    漫谈边缘计算(二):各怀心事的玩家
    漫谈边缘计算(一):边缘计算是大势所趋
    从小小后视镜看物联网的生态(下)
    机器学习笔记(四)---- 逻辑回归的多分类
    机器学习笔记(三)---- 逻辑回归(二分类)
    机器学习笔记(二)---- 线性回归
  • 原文地址:https://www.cnblogs.com/cjyyb/p/9268527.html
Copyright © 2011-2022 走看看