自然数幂和

自然数幂和

一个经典的问题是求自然数幂和：
\(\sum_{i=1}^{n}i^k=1^k+2^k+...+n^k\)
根据伯努利的公式：
\(\sum_{i=1}^{n}i^k=\frac{1}{k+1}\sum_{i=1}^{k+1}C_{k+1}^{i}B_{k+1-i}(n+1)^i\)条件：\(B_1=-\frac{1}{2}\)
知道前k+1的伯努利数之后我们就可以\(O(k)\)地求出自然数幂和
但是伯努利数并不好求。

1.根据\(B_n=-\frac{1}{n+1}\sum_{i=0}^{n-1}C_{n+1}^{i}B_{i}\)
我们可以\(O(k^2)\)的求出
2.根据母函数

\[\sum_{i=0}^{\infty}\frac{B_{i}}{i!}x^i=\frac{x}{e^x-1}=\frac{1}{\sum_{i=0}^{\infty}\frac{x^i}{(i+1)!}} \]

求前k项，只需求\({\sum_{i=0}^{\infty}\frac{x^i}{(i+1)!}}\ mod\ x^{k+1}的逆元\)

3.设\(S(n,k)=\sum_{i=1}^{n}i^k\)
\(S(n,k)=\sum_{i=0}^{k+1}S(i,k)(-1)^{k+1-i}\frac{\prod_{j=0}^{k+1}n-j}{(n-i)i!(k+1-i)!}\)
可以O(k)的求出

下面有2、3的代码

代码

2：
(http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1258)
由于求的是模1e9+7，不能用NTT,只能用取模的FFT来算，虽然我感觉常数巨大，但相比T*k可以忽略，跑得反而是最快的，比3还快，大概是因为询问的时候常数很小。

//#include <bits/stdc++.h>
#include <stdio.h>
#include <iostream>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <algorithm>
#include <queue>
#include <vector>
#include <set>
#include <map>
#include <stack>
#include <bitset>
#include <string>
#include <time.h>
using namespace std;
long double esp=1e-11;
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#define fi first
#define se second
#define all(a) (a).begin(),(a).end()
#define cle(a) while(!a.empty())a.pop()
#define mem(p,c) memset(p,c,sizeof(p))
#define mp(A, B) make_pair(A, B)
#define pb push_back
#define lson l , m , rt << 1
#define rson m + 1 , r , rt << 1 | 1
typedef long long int LL;
const long double PI = acos((long double)-1);
const LL INF=0x3f3f3f3fll;
const int MOD =1000000007ll;
const int maxn=140100;
struct complex
{
    double r,i;
    complex(double _r = 0.0,double _i = 0.0){r = _r; i = _i;}
    complex operator +(const complex &b){return complex(r+b.r,i+b.i);}
    complex operator -(const complex &b){return complex(r-b.r,i-b.i);}
    complex operator *(const complex &b){return complex(r*b.r-i*b.i,r*b.i+i*b.r);}
};
complex conj(complex a){return complex(a.r,-a.i);}
complex w[maxn];
int bitrev[maxn];
inline void fft_prepare(int len)
{
	int L=__builtin_ctz(len);
	for(int i=0;i<len;i++) bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (L - 1));
	for(int i=0;i<len;i++) w[i] = complex(cos(2 * PI * i / len), sin(2 * PI * i / len));
}
void FFT(complex *a, const int &n)
{
	for(int i=0;i<n;i++) if (i < bitrev[i]) swap(a[i], a[bitrev[i]]);
	for (int i = 2, lyc = n >> 1; i <= n; i <<= 1, lyc >>= 1)
		for (int j = 0; j < n; j += i)
		{
			complex *l = a + j, *r = a + j + (i >> 1), *p = w;
			for(int k=0;k<(i >> 1);k++)
			{
				complex tmp = *r * *p;
				*r = *l - tmp, *l = *l + tmp;
				++l, ++r, p += lyc;
			}
		}
}
int callen(LL len1,LL len2){
    int len=1;
    while(len < (len1<<1) || len < (len2<<1))len<<=1;
    return len;
}
LL fftans[maxn];
complex A[maxn],B[maxn],dft[4][maxn],dt[4];
LL td[4];
void fft(LL* y1,int len1,LL * y2,int len2,LL mod){
    int len=callen(len1,len2);
    fft_prepare(len);
    for(int x=0;x<len1;x++)A[x]=complex(y1[x]&32767,y1[x]>>15);
    for(int x=len1;x<len;x++)A[x]=complex(0,0);
    for(int x=0;x<len2;x++)B[x]=complex(y2[x]&32767,y2[x]>>15);
    for(int x=len2;x<len;x++)B[x]=complex(0,0);
    FFT(A,len);FFT(B,len);
    int j;
    for(int x=0;x<len;x++)
    {
        j=(len-x)&(len-1);
        dt[0]=(A[x]+conj(A[j]))*complex(0.5,0);
        dt[1]=(A[x]-conj(A[j]))*complex(0,-0.5);
        dt[2]=(B[x]+conj(B[j]))*complex(0.5,0);
        dt[3]=(B[x]-conj(B[j]))*complex(0,-0.5);
        dft[0][j]=dt[0]*dt[2];
        dft[1][j]=dt[0]*dt[3];
        dft[2][j]=dt[1]*dt[2];
        dft[3][j]=dt[1]*dt[3];
    }
    for(int x=0;x<len;x++)
    {
        A[x]=dft[0][x]+dft[1][x]*complex(0,1);
        B[x]=dft[2][x]+dft[3][x]*complex(0,1);
    }
    FFT(A,len);FFT(B,len);
    for(int x=0;x<len;x++)
    {
        td[0]=(LL)(A[x].r/len+0.5)%mod;
        td[1]=(LL)(A[x].i/len+0.5)%mod;
        td[2]=(LL)(B[x].r/len+0.5)%mod;
        td[3]=(LL)(B[x].i/len+0.5)%mod;
        fftans[x]=(td[0]+((LL)(td[1]+td[2])<<15)+((LL)td[3]<<30))%mod;
    }
}
LL Inv(LL a,LL p) //0<a<p,a%=p先
{
    if(a==0)return -1;
    if(a==1)return 1;
    return Inv(p%a,p)*(p-p/a)%p;
}
void getinv(LL C[],LL D[],int t)
{
	if(t==1){D[0]=Inv(C[0],MOD);return;}
	getinv(C,D,(t+1)>>1);
	fft(C,t,D,t,MOD);
	fft(fftans,t,D,t,MOD);
	for(int x=0;x<t;x++){D[x]=(2ll*D[x]-fftans[x])%MOD;if(D[x]<0)D[x]+=MOD;}
	
}
LL o[maxn],Ans[maxn],Ber[maxn];
LL inv[maxn],pin_inv[maxn],pin[maxn];
void initCab(int N){
	inv[0]=pin[0]=pin_inv[0]=inv[1]=pin[1]=pin_inv[1]=1;
	for(int x=2;x<=N;x++){
		inv[x]=(LL)(MOD-MOD/x)*inv[MOD%x]%MOD;
		pin[x]=(LL)pin[x-1]*x%MOD;
		pin_inv[x]=(LL)pin_inv[x-1]*inv[x]%MOD;}}
LL Cab(int n,int m){return (LL)pin[n]*pin_inv[m]%MOD*pin_inv[n-m]%MOD;}
int main()
{
    //freopen("in.txt", "r", stdin);
    //freopen("inlay.in", "r", stdin);
    //freopen("out.txt", "w", stdout);
    //::iterator iter;                  %I64d
    //for(int x=1;x<=n;x++)
    //for(int y=1;y<=n;y++)
    //scanf("%d",&a);
    //printf("%d\n",ans);
    int sz=50000;
    initCab(sz+10);			
    o[0]=1;
    for(int x=1;x<=sz;x++)
		o[x]=pin_inv[x+1];
	getinv(o,Ans,sz+1);
	for(int x=0;x<=sz;x++)
		Ber[x]=(LL)Ans[x]*pin[x]%MOD;
	int T;			
	scanf("%d",&T);
	while(T--)
	{
		LL n,k;
		scanf("%lld%lld",&n,&k);
		n%=MOD;
		LL ans=0,tem=1;
		for(int x=1;x<=k+1;x++)
		{
			tem=tem*(n+1)%MOD;
			ans=(ans+1ll*Cab(k+1,x)*Ber[k+1-x]%MOD*tem)%MOD;
		}
		ans=ans*inv[k+1]%MOD;
		printf("%lld\n",ans);
	}
    return 0;
}

LL mu(LL a,LL b)		//a负数注意
{
    LL ans=1;
    while(b)
    {
        if((b & 1)&&(ans*=a)>=MOD)
            ans%=MOD;
        b >>= 1;
        if((a*=a)>=MOD)
            a%=MOD;
    }
    return ans;
}
LL pin[maxn],prime[maxn],bo[maxn],cnt,S[maxn],L[maxn],R[maxn];
void initsp(int sz)
{
	pin[0]=1;
	for(int x=1;x<=sz;x++)
		pin[x]=pin[x-1]*mu(x,MOD-2)%MOD;
	for(int x=2;x<=sz;x++){
		if(!bo[x])prime[++cnt]=x;
		for(int y=1;y<=cnt&&x*prime[y]<=sz;y++)
		{
			bo[x*prime[y]]=prime[y];
			if(x%prime[y])break;
		}
	}
}
LL sumpow(LL n,LL k){         //模意义下的求和i^k	O(k) 0^0=0
	LL ans=0,sum=0;
	for(int x=1;x<=k+1;x++)if(!bo[x])S[x]=mu(x,k);else S[x]=S[x/bo[x]]*S[bo[x]]%MOD;
	for(int x=1;x<=k+1;x++)S[x]=(S[x]+S[x-1])%MOD;
	for(int x=0;x<=k+1;x++)L[x]=R[x]=n-x;
	for(int x=1;x<=k+1;x++)L[x]=L[x-1]*L[x]%MOD;
	for(int x=k;x>=0;x--)R[x]=R[x+1]*R[x]%MOD;
	for(int x=0;x<=k+1;x++){
		sum=S[x];
		if(x>0)sum=sum*pin[x]%MOD*L[x-1]%MOD;
		if(x<k+1)sum=sum*pin[k+1-x]%MOD*((k+1-x)%2?MOD-1:1)%MOD*R[x+1]%MOD;
		ans=(ans+sum)%MOD;
	}
	return (ans+MOD)%MOD;
}