自然数幂和
一个经典的问题是求自然数幂和：
(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
",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
",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;
}