题目概述

求 $\sum_{i=1}^{n}i^K,K\le 50000$ 。

解题报告

同51Nod1228，只不过 $K\le 50000$ 所以不能 $O(K^2)$ 处理伯努利数，于是我们根据 $\hat{B}(x)={x\over e^x-1}$ 和多项式求逆来预处理，把 $x$ 除下去，其实就是 $e^x$ 往左移动了一位，只需要改变下下标：

$$ \hat{B}(x)=(\sum_{i\ge0}{x^i\over(i+1)!})^{-1} $$

对后面那个多项式进行求逆，得到 $\hat{B}(x)$ ，那么 $B(x)_i=\hat{B}(x)_ii!$ 。

因为这道题不是NTT模数所以需要三模数NTT或者拆系数FFT……

示例程序

#include<cmath>
#include<cstdio>
#include<algorithm>
#define fr first
#define sc second
#define mp make_pair
using namespace std;
typedef long long LL;typedef long double DB;typedef pair<DB,DB> CN;
const int maxk=50000,maxt=1<<17,MOD=1e9+7;const DB pi=acos(-1);

int te,K;LL n;int INV[maxt+5],pre[maxt+5],fac[maxt+5],A[maxt+5],Be[maxt+5],tem[maxt+5],pw[maxk+5];
int rev[maxt+5];DB Sin[maxt+5],Cos[maxt+5];

inline void Pre(int n){
    for (int i=1;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)*(n>>1));
    for (int i=1;i<n;i<<=1) Sin[i]=sin(pi/i),Cos[i]=cos(pi/i);
}
inline CN operator + (const CN &a,const CN &b) {return mp(a.fr+b.fr,a.sc+b.sc);}
inline CN operator - (const CN &a,const CN &b) {return mp(a.fr-b.fr,a.sc-b.sc);}
inline CN operator * (const CN &a,const CN &b) {return mp(a.fr*b.fr-a.sc*b.sc,a.fr*b.sc+a.sc*b.fr);}
inline void FFT(CN *a,int n,int f){
    for (int i=0;i<n;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int k=1;k<n;k<<=1){
        CN wn=mp(Cos[k],f>0?Sin[k]:-Sin[k]),w=mp(1,0),x,y;
        for (int i=0;i<n;i+=k<<1,w=mp(1,0))
            for (int j=0;j<k;j++,w=w*wn)
                x=a[i+j],y=a[i+j+k]*w,a[i+j]=x+y,a[i+j+k]=x-y;
    }if (f<0) for (int i=0;i<n;i++) a[i].fr/=n;
}
inline void AMOD(int &x,int tem) {if ((x+=tem)>=MOD) x-=MOD;}
inline void MTT(int *a,int *b,int *ans,int n){
    static CN A[maxt+5],B[maxt+5],C[maxt+5],D[maxt+5],E[maxt+5],F[maxt+5],G[maxt+5],H[maxt+5];
    for (int i=0;i<n;i++) A[i]=mp(a[i]>>15,0),B[i]=mp(a[i]&32767,0),A[i+n]=B[i+n]=mp(0,0);
    for (int i=0;i<n;i++) C[i]=mp(b[i]>>15,0),D[i]=mp(b[i]&32767,0),C[i+n]=D[i+n]=mp(0,0);
    int m=n<<1;Pre(m);FFT(A,m,1);FFT(B,m,1);FFT(C,m,1);FFT(D,m,1);
    for (int i=0;i<m;i++) E[i]=A[i]*C[i],F[i]=A[i]*D[i],G[i]=B[i]*C[i],H[i]=B[i]*D[i];
    FFT(E,m,-1);FFT(F,m,-1);FFT(G,m,-1);FFT(H,m,-1);for (int i=0;i<n;i++) ans[i]=0;
    for (int i=0;i<n;i++){
        AMOD(ans[i],((LL)(E[i].fr+0.5)%MOD<<30)%MOD);AMOD(ans[i],((LL)(F[i].fr+0.5)%MOD<<15)%MOD);
        AMOD(ans[i],((LL)(G[i].fr+0.5)%MOD<<15)%MOD);AMOD(ans[i],(LL)(H[i].fr+0.5)%MOD);
    }
}
inline int Pow(int w,int b) {int s;for (s=1;b;b>>=1,w=(LL)w*w%MOD) if (b&1) s=(LL)s*w%MOD;return s;}
void Inv(int *A,int *B,int n){
    if (n==1) {B[0]=Pow(A[0],MOD-2);return;}Inv(A,B,n>>1);MTT(B,B,tem,n);MTT(A,tem,tem,n);
    for (int i=0;i<n;i++) AMOD(B[i],B[i]),AMOD(B[i],MOD-tem[i]);
}
inline void Make(int n){
    int t;for (t=1;t<=n;t<<=1);pre[0]=INV[1]=1;fac[0]=1;
    for (int i=2;i<=t;i++) INV[i]=(LL)(MOD-MOD/i)*INV[MOD%i]%MOD;
    for (int i=1;i<=t;i++) pre[i]=(LL)pre[i-1]*INV[i]%MOD,fac[i]=(LL)fac[i-1]*i%MOD;
    for (int i=0;i<t;i++) A[i]=pre[i+1];Inv(A,Be,t);for (int i=0;i<=n;i++) Be[i]=(LL)Be[i]*fac[i]%MOD;
}
#define C(x,y) ((x)<(y)?0:((LL)fac[x]*pre[y]%MOD*pre[(x)-(y)]%MOD))
int main(){
    freopen("program.in","r",stdin);freopen("program.out","w",stdout);
    for (Make(maxk+1),scanf("%d",&te);te;te--){
        scanf("%lld%d",&n,&K);n++;n%=MOD;int ans=0;
        pw[0]=1;for (int i=1;i<=K+1;i++) pw[i]=(LL)pw[i-1]*n%MOD;
        for (int i=0;i<=K;i++) AMOD(ans,(LL)C(K+1,i)*Be[i]%MOD*pw[K+1-i]%MOD);
        ans=(LL)ans*INV[K+1]%MOD;printf("%d\n",ans);
    }return 0;
}