51NOD 1258 自然数幂和

题目链接

description

$T$ 次询问，每次给定 $n,k$，求 $\sum _{i=1}^n{i}^k$ 模 1e9+7.

$n\le {10}^{18},k\le 5\times {10}^4$

solution

可以拉插用什么多项式

考虑将 $n$ 带入 $f(x)=\sum _{i=1}^x{i}^k$。可以证明，$f(x)$ 是一个 $k+1$ 次多项式。

一种感性的理解方式是 $i^k$ 积分后是个 $k+1$ 次函数

于是我们可以通过拉格朗日插值法构造出这个多项式在 $n$ 处的取值。

我们容易线性筛求出 $f(x)$ 在 $1,2,\dots ,k,k+1,k+2$ 处的点值，根据 $n+1$ 个点确定一个 $n$ 次多项式，我们可以通过这些点值构造出这个多项式。

设 $t=k+2$，$s_p=\sum _{i=1}^p{i}^k$，对于 $1\le p\le t,p\in \mathbb{Z}$，构造函数 $g_p(x)=s_p\prod _{i=1,i\ne p}^t{\displaystyle \frac{(x-i)}{(p-i)}}$。不难验证，$\mathrm{\forall }x\in [1,t]\cap \mathbb{Z},x\ne p$ 有 $g_p(x)=0$ 且 $g_p(p)=s_p$。

于是可以得出 $f(x)=\sum _{p=1}^t{g}_p(x)=\sum _{p=1}^t{s}_p\prod _{i=1,i\ne p}^t{\displaystyle \frac{(x-i)}{(p-i)}}$。

于是 $f(n)=\sum _{p=1}^t{s}_p\prod _{i=1,i\ne p}^t{\displaystyle \frac{(n-i)}{(p-i)}}$

线性筛出 $s_p$，预处理阶乘以及 $n-i$ 的前缀后缀积，即可做到严格线性。

时间复杂度 $O(Tk)$

code

#include<bits/stdc++.h>
 
using namespace std;
 
const int N=(1<<18)+10,mod=1e9+7;
 
int f[N],n,k,m;
bitset<N> st;
 
int ksm(int a,int b){
  int ret=1;
  while(b){
    if(b&1) ret=ret*1ll*a%mod;
    a=a*1ll*a%mod;
    b>>=1;
  }
  return ret;
}
 
vector<int> primes;
 
int pre[N],sfx[N],jc[N];
 
void solve(){
  long long nn;
  cin>>nn;
  nn%=mod;
  cin>>k;
  n=nn;
 
  primes.clear();
  st.reset();
  m=k+2;
  f[1]=1;
  for(int i=2; i<=m; i++){
    if(!st[i]){
      f[i]=ksm(i,k);
      primes.emplace_back(i);
    }
    for(auto j:primes){
      if(i*j>m) break;
      st[i*j]=1;
      f[i*j]=f[i]*1ll*f[j]%mod;
      if(i%j==0){
        break;
      }
    }
  }
  for(int i=2; i<=m; i++){
    f[i]=(f[i-1]+f[i])%mod;
  }
 
  jc[0]=pre[0]=sfx[m+1]=1;
  for(int i=1; i<=m; i++){
    jc[i]=jc[i-1]*1ll*i%mod;
    pre[i]=pre[i-1]*1ll*(n-i)%mod;
  }
 
  for(int i=m; i; i--){
    sfx[i]=sfx[i+1]*1ll*(n-i)%mod;
    jc[i]=i==m?ksm(jc[i],mod-2):jc[i+1]*1ll*(i+1)%mod;
  }
 
  int ans=0;
  for(int i=1; i<=m; i++){
    int t=sfx[i+1]*1ll*pre[i-1]%mod*jc[i-1]%mod*jc[m-i]%mod;
    if((m-i)&1) t=mod-t;
    ans=(ans+t*1ll*f[i]%mod)%mod;
  }
 
  cout<<ans<<'\n';
}
 
int main(){
 
  int T;
  cin>>T;
  while(T--) solve();
 
  return 0;
}