51Nod 1258 序列求和 V4（组合数学+FFT）

Description

$T(n) = n^k，S(n) = T(1) + T(2) + ...... T(n)$。给出$n$和$k$，求$S(n)$。

例如$k = 2，n = 5，S(n) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 = 55$。

由于结果很大，输出$S(n)\ Mod\ 1000000007$的结果即可。

Input

第$1$行：一个数$T$，表示后面用作输入测试的数的数量。
第$2 - T + 1$行：每行$2$个数，$n, k$中间用空格分割。

$(1\le T\le 500,1\le n\le 10^{18},1\le k\le 5\cdot 10^4)$

Output

共$T$行，对应$S(n)\ Mod\ 1000000007$的结果。

Sample Input

3
5 3
4 2
4 1

Sample Output

225
30
10

Solution

求自然数幂和，即$S_k(n)=\sum\limits_{i=1}^ni^k$，则有$S_k(n)=\frac{1}{k+1}\sum\limits_{i=0}^{k+1}C_{k+1}^iB_{k+1-i}(n+1)^i$，伯努利数$O(k^2)$预处理显然不行，但由于$\sum\limits_{i=0}^nC_{n+1}^iB_i=(n+1)!\sum\limits_{i=0}^n\frac{B_i}{i!}\frac{1}{(n+1-i)!}=0$，而$B_0=1$，故多项式$f(x)=\sum\limits \frac{B_i}{i!}x^i$即为多项式$g(x)=\sum\limits \frac{(n+1-i)!}x^i$的逆，用多项式求逆即可$O(klog^2k)$预处理出$B_0,...,B_{k+1}$，单次查询时间复杂度$O(k)$

Code

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
#define maxfft 131072+5
#define mod 1000000007
int add(int x,int y)
{
    x+=y;
    if(x>=mod)x-=mod;
    return x;
}
int mul(int x,int y)
{
    ll z=1ll*x*y;
    return z-z/mod*mod;
}
int Pow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=mul(ans,a);
        a=mul(a,a);
        b>>=1;
    }
    return ans;
}
const double pi=acos(-1.0);
struct cp
{
    double a,b;
    cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}
    cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}
    cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}
    cp operator *(const double &o)const {return (cp){a*o,b*o};}
    cp operator !() const{return (cp){a,-b};}
}w[maxfft];
int pos[maxfft];
void fft_init(int len)
{
    int j=0;
    while((1<<j)<len)j++;
    j--;
    for(int i=0;i<len;i++)
        pos[i]=pos[i>>1]>>1|((i&1)<<j);
}
void fft(cp *x,int len,int sta)
{
    for(int i=0;i<len;i++)
        if(i<pos[i])swap(x[i],x[pos[i]]);
    w[0]=(cp){1,0};
    for(int i=2;i<=len;i<<=1)
    {
        cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};
        for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];
        for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;
        for(int j=0;j<len;j+=i)
        {
            cp *a=x+j,*b=a+(i>>1);
            for(int l=0;l<i>>1;l++)
            {
                cp o=b[l]*w[l];
                b[l]=a[l]-o;
                a[l]=a[l]+o;
            }
        }
    }
    if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;
}
cp x[maxfft],y[maxfft],z[maxfft];
void FFT(int *a,int *b,int n,int m,int *c)
{
    int len=1;
    while(len<n+m)len<<=1;
    fft_init(len);
    for(int i=0;i<len;i++)
    {
        int aa=i<n?a[i]:0,bb=i<m?b[i]:0;
        x[i]=(cp){(aa>>15),(aa&32767)},y[i]=(cp){(bb>>15),(bb&32767)};
    }
    fft(x,len,1),fft(y,len,1);
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=((x[i]+!x[j])*(y[i]-!y[j])+(x[i]-!x[j])*(y[i]+!y[j]))*(cp){0,-0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod;
        ta=(ta<<15)%mod;
        c[i]=ta;
    }
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=(x[i]-!x[j])*(y[i]-!y[j])*(cp){-0.25,0}+(x[i]+!x[j])*(y[i]+!y[j])*(cp){0,0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod,tb=(ll)(z[i].b+0.5)%mod;
        ta=(ta+(tb<<30))%mod;
        c[i]=(c[i]+ta)%mod;
    }
}
int temp1[maxfft];
void Poly_Inv(int *poly,int n,int *ans)
{
    ans[0]=Pow(poly[0],mod-2);
    for(int i=2;i<=n;i<<=1)
    {
        FFT(poly,ans,i,i/2,temp1);
        FFT(ans,temp1+i/2,i/2,i/2,temp1);
        for(int j=0;j<i/2;j++)ans[j+i/2]=temp1[j]==0?0:mod-temp1[j];
    }
}
int inv[maxfft],fact[maxfft],B[maxfft];
void init(int n=50001)
{
    fact[0]=1;
    for(int i=1;i<=n;i++)fact[i]=mul(i,fact[i-1]);
    inv[1]=1;
    for(int i=2;i<=n;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
    inv[0]=1;
    for(int i=1;i<=n;i++)inv[i]=mul(inv[i-1],inv[i]);
    int len=1;
    while(len<n)len<<=1;
    Poly_Inv(inv+1,len,B);
    for(int i=0;i<=n;i++)B[i]=mul(B[i],fact[i]);
}
int C(int n,int m)
{
    if(m<0||m>n)return 0;
    return mul(fact[n],mul(inv[m],inv[n-m]));
}
int main()
{
    init();
    int T,k;
    ll n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%d",&n,&k);
        int ans=0;
        for(int i=1;i<=k+1;i++)ans=add(ans,mul(mul(C(k+1,i),B[k+1-i]),Pow((n+1)%mod,i)));
        ans=mul(ans,Pow(k+1,mod-2));
        printf("%d\n",ans);
    }
}