• 伯努利数学习笔记


    1.定义式

    定义伯努利数列(B_n)满足:

    [B_0=1,sum_{i=0}^n{n+1choose i}B_i=0(n>0) ]

    2.递推式

    可以发现定义式里面包含了(B_n)这一项,于是把(B_n)提出来:

    [-{n+1choose n}B_n=sum_{i=0}^{n-1}{n+1choose i}B_i \-(n+1)B_n=sum_{i=0}^{n-1}{n+1choose i}B_i \B_n=-frac{1}{n+1}sum_{i=0}^{n-1}{n+1choose i}B_i ]

    直接用定义式求是(O(n^2))的复杂度

    3.生成函数

    把定义式的循环上界减一,得:

    [sum_{i=0}^{n-1}{nchoose i}B_i=0 ]

    注意到组合数上标变成了(n),再加个(B_n)

    [sum_{i=0}^{n-1}{nchoose i}B_i+B_n=B_n \sum_{i=0}^{n-1}{nchoose i}B_i+{nchoose n}B_n=B_n \sum_{i=0}^{n}{nchoose i}B_i=B_n ]

    组合数很烦,把它拆开来:

    [sum_{i=0}^{n}frac{n!}{i!(n-i)!}B_i=B_n \sum_{i=0}^nfrac{B_i}{i!}frac{1}{(n-i)!}=frac{B_n}{n!} ]

    两边都写成生成函数的形式:

    [sum_{n=0}(sum_{i=0}^nfrac{B_i}{i!}frac{1}{(n-i)!})x^n =sum_{n=0}(frac{B_n}{n!})x^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)项:

    [[x^n]B(x)e^x=sum_{i=0}^nfrac{B_i}{i!}frac{1}{(n-i)!} \ [x^0]B(x)e^x=B_0 \ [x^1]B(x)e^x=B_0+B_1 ]

    (B(x))的第一项只有(B_1x),少了个(B_0x),加上去即可:

    [B(x)e^x=B(x)+B_0x=B(x)+x ]

    可以解出(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))的指数型生成函数为

    [S_m(x)=sum_{k=0}frac{sum_{i=1}^mi^k}{k!}x^k ]

    交换求和符号得到:

    [S_m(x)=sum_{i=1}^msum_{k=0}frac{i^k}{k!}x^k \=sum_{i=1}^me^{ix}=frac{e^{(m+1)x}-e^x}{e^x-1} ]

    注意到分母和(B(x))一样,把(B(x))带进去:

    [S_m(x)=B(x)frac{e^{(m+1)x}-e^x}{x}=B(x)frac{e^x}{x}(e^{mx}-1) ]

    (e^{mx})展开得到:

    [S_m(x)=B(x)frac{e^x}{x}(sum_{i=0}frac{(mx)^i}{i!}-1) \=B(x)frac{e^x}{x}sum_{i=1}frac{(mx)^i}{i!} \=B(x)e^xsum_{i=1}frac{m^ix^{i-1}}{i!} \=B(x)e^xsum_{i=0}frac{m^{i+1}x^i}{(i+1)!} ]

    前面这个(B(x)e^x)不好处理,但注意到(B(x)e^x=B(x)+x),于是引入新的(B'(x)=B(x)+x),把(B'(x))展开:

    [S_m(x)=(sum_{i=0}frac{B'_i}{i!}x^i)(sum_{i=0}frac{m^{i+1}}{(i+1)!}x^i) \=sum_{n=0}(sum_{i=0}^nfrac{B'_i}{i!}frac{m^{n-i+1}}{(n-i+1)!})x^n \=sum_{n=0}(sum_{i=0}^n{n+1choose i}B'_im^{n-i+1})frac{x^n}{(n+1)!} \=sum_{n=0}(frac{1}{n+1}sum_{i=0}^n{n+1choose i}B'_im^{n-i+1})frac{x^n}{n!} ]

    注意到(S_m(x)=sum_{k=0}S(m,k)frac{x^k}{k!}),所以:

    [S(m,k)=frac{1}{k+1}sum_{i=0}^k{k+1choose i}B'_im^{k-i+1} ]

    这个就可以称为自然数幂和的通项公式
    注意到通项公式可以写成卷积的形式:

    [S(m,k)=frac{1}{k+1}sum_{i=0}^k{k+1choose i}B'_im^{k-i+1} \=k!sum_{i=0}^kfrac{B'_i}{i!}frac{m^{k-i+1}}{(k-i+1)!} ]

    因为(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),这样答案的式子就是:

    [a_0+sum_{k=0}^nS_k(x)a_k ]

    只需要关注前面的求和式,用伯努利数将其展开得到:

    [sum_{k=0}^n(k!sum_{i=0}^kfrac{B'_i}{i!}frac{x^{k-i+1}}{(k-i+1)!})a_k \=sum_{k=0}^na_kk!sum_{i=0}^kfrac{B'_i}{i!}frac{x^{k-i+1}}{(k-i+1)!} \=sum_{k=0}^na_kk!sum_{i=0}^kfrac{B'_{k-i}}{(k-i)!}frac{x^{i+1}}{(i+1)!} \=sum_{k=0}^na_kk!sum_{i=1}^{k+1}frac{B'_{k-i+1}}{(k-i+1)!}frac{x^i}{i!} \=sum_{i=1}^{n+1}frac{x^i}{i!}sum_{k=i-1}^na_kk!frac{B'_{k-i+1}}{(k-i+1)!} ]

    那么设(f_i=a_ii!)(g_i=frac{B'_i}{i!}),引入(g'_i=g_{n-i}),那么上面的式子可以接着写成:

    [sum_{i=1}^{n+1}frac{x^i}{i!}sum_{k=i-1}^nf_kg'_{n-k+i-1} ]

    注意到(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;
    }
    
  • 相关阅读:
    如何编译树莓派内核
    代码导出Reporting Services报表文件
    Bit-Coin收入的一分钱
    如何在树莓派上运行雷神之锤III
    新树莓派入手
    如何通过PowerShell在Visual Studio的Post-build中预热SharePoint站点
    每日一题20201218(389. 找不同)
    每日一题20201217(714. 买卖股票的最佳时机含手续费)
    每日一题20201216(290. 单词规律)
    每日一题20201215(738. 单调递增的数字)
  • 原文地址:https://www.cnblogs.com/akura/p/12577575.html
Copyright © 2020-2023  润新知