CFP622F The Sum of the k-th Powers 学习笔记

CFP622F The Sum of the k-th Powers 学习笔记

Luogu Link

题意简述

求 \(\sum_{i=1}^n i^k\)。答案对 \(10^9+7\) 取模。

\(n\le 10^9\)，\(k\le 10^6\)。

做法解析

题目即问我们“自然数幂之和”，它可以记作 \(S_{k}(n)=\sum_{i=1}^n i^k\)。

怎么求呢？首先我们发现，\(n\) 在 \(10^9\) 范围内，这说明我们的复杂度里不能出现任何不低于一次的 \(n\) 项，怎么办呢？

于是我们考虑将这个式子看作一个 \(n\) 作为自变量的函数去思考。正好，这个多项式的次数是 \(k+1\) 次的。所以它可以通过拉插 \(k+2\) 个点值来唯一确定。

为什么它是 \(k+1\) 次的多项式？

对于一个数列 \(\{a_n\}\)，我们作其相邻元素的差分得到一个新的数列 \(\{b_n\}\)，我们称 \(\{b_n\}\) 为 \(\{a_n\}\) 的一阶差分数列。\(\{a_n\}\) 的二阶差分就是 \(\{b_n\}\) 的一阶差分。依此类推。如果数列 \(\{a_n\}\) 的 \(p\) 阶差数列是个非 \(0\) 常数数列，那么称它为 \(p\) 阶等差数列。

又这里有个定理：数列 \(\{a_n\}\) 是 \(p\) 阶等差数列的充要条件是数列的通项 \(a_n\) 是关于 \(n\) 的一个 \(p\) 次多项式。因为 \(\{a_n\}\) 每做一次差分，最高次项都会因为 \(f_i((x+1)^k-x^k)=f_i(x^k+\cdots+1-x^k)\) 被消掉。

所以你发现，当我们给 \(a_n=\sum_{i=1}^n x^k\) 做差分的时候，刚好有 \(b_n=n^k\)，而这是关于 \(n\) 的一个 \(k\) 次多项式。所以原柿是关于 \(n\) 的一个 \(k+1\) 次多项式。

看懂了吗？

还是不理解？你这么看：原柿子还可以写成 \(\sum_{i=1}^n (n-i)^k\)，显然拆开后会出现 \(n\) 次 \(n^k\)，这不就是 \(n^{k+1}\)？。

拉插一般而言是 \(O(k^2)\) 的。不过这道题所有点值的自变量是连续等差 \(1\) 的，所以有 \(O(k)\) 方法可以计算。详见LGP4718学习笔记。就这样。

代码实现

#include <bits/stdc++.h>
using namespace std;
using namespace obasic;
using namespace omodint;
using mint=m107;
const int MaxK=1e6+5;
mint facr[MaxK],finv[MaxK];
mint pprd[MaxK],sprd[MaxK];
namespace omathe{
    void premwork(int n){
        facr[0]=finv[0]=1;for(int i=1;i<=n;i++)facr[i]=facr[i-1]*i;
        finv[n]=1/facr[n];for(int i=n-1;i>=1;i--)finv[i]=finv[i+1]*(i+1);
    }
    mint fpow(mint a,int b){
        mint res=1;
        for(;b;a*=a,b>>=1)if(b&1)res*=a;
        return res;
    }
    mint laiitp(mint a[],int n,mint k){
        mint res=0;
        pprd[0]=k;for(int j=1;j<=n;j++)pprd[j]=pprd[j-1]*(k-j);
        sprd[n]=k-n;for(int j=n-1;j>=0;j--)sprd[j]=sprd[j+1]*(k-j);
        for(int i=0;i<=n;i++)res+=a[i]*(i>0?pprd[i-1]:1)*(i<n?sprd[i+1]:1)*finv[i]*finv[n-i]*((n-i)&1?-1:1);
        return res;
    }
};
using namespace omathe;
int N,K;mint A[MaxK];
int main(){
    readis(N,K);premwork(K+1);
    for(int i=1;i<=K+1;i++)A[i]=A[i-1]+fpow(i,K);
    writi(miti(laiitp(A,K+1,N)));
    return 0;
}