Description
T(n)=nk,S(n)=T(1)+T(2)+......T(n) T ( n ) = n k , S ( n ) = T ( 1 ) + T ( 2 ) + . . . . . . T ( n ) 。给出 n n 和,求 S(n) S ( n ) 。
例如 k=2,n=5,S(n)=12+22+32+42+52=55 k = 2 , n = 5 , S ( n ) = 1 2 + 2 2 + 3 2 + 4 2 + 5 2 = 55 。
由于结果很大,输出 S(n) Mod 1000000007 S ( n ) M o d 1000000007 的结果即可。
Input
第
1
1
行:一个数,表示后面用作输入测试的数的数量。
第
2−T+1
2
−
T
+
1
行:每行
2
2
个数,中间用空格分割。
(1≤T≤500,1≤n≤1018,1≤k≤5⋅104) ( 1 ≤ T ≤ 500 , 1 ≤ n ≤ 10 18 , 1 ≤ k ≤ 5 ⋅ 10 4 )
Output
共 T T 行,对应的结果。
Sample Input
3
5 3
4 2
4 1
Sample Output
225
30
10
Solution
求自然数幂和,即 Sk(n)=∑i=1nik S k ( n ) = ∑ i = 1 n i k ,则有 Sk(n)=1k+1∑i=0k+1Cik+1Bk+1−i(n+1)i S k ( n ) = 1 k + 1 ∑ i = 0 k + 1 C k + 1 i B k + 1 − i ( n + 1 ) i ,伯努利数 O(k2) O ( k 2 ) 预处理显然不行,但由于 ∑i=0nCin+1Bi=(n+1)!∑i=0nBii!1(n+1−i)!=0 ∑ i = 0 n C n + 1 i B i = ( n + 1 ) ! ∑ i = 0 n B i i ! 1 ( n + 1 − i ) ! = 0 ,而 B0=1 B 0 = 1 ,故多项式 f(x)=∑Bii!xi f ( x ) = ∑ B i i ! x i 即为多项式 g(x)=∑(n+1−i)!xi g ( x ) = ∑ ( n + 1 − i ) ! x i 的逆,用多项式求逆即可 O(klog2k) O ( k l o g 2 k ) 预处理出 B0,...,Bk+1 B 0 , . . . , B k + 1 ,单次查询时间复杂度 O(k) 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);
}
}