[题解] CF622F The Sum of the k-th Powers

Posted 2021-03-13 wxq1229

CF622F The Sum of the k-th Powers

题意：给(n)和(k)，让你求(sumlimits_{i = 1} ^ n i^k mod 10^9 + 7)。((1 le n le 10^9,0 le k le 10^6))

---

好，我们先不看题，来补一些数学。

考虑这样一个序列
[ h_0,h_1,dots,h_n,dots ]
我们定义它的一个差分序列（一阶）
[ Delta h_0, Delta h_1, dots , Delta h_n, dots ]

满足(Delta h_n = h_{n+1} - h_n (n ge 0))，换句话说（一阶）差分序列就是原序列相邻两项的差。

同样的，我们还能在原序列的一阶差分序列上再做一次差分，得到二阶差分序列：
[ Delta(Delta h_0), Delta(Delta h_1),dots , Delta(Delta h_n), dots ]
可以把前面的(Delta)写成(Delta ^2)，就是(Delta^2 h_0, Delta^2 h_1,dots , Delta^2 h_n, dots)

根据定义，有(Delta^2 h_n = Delta(Delta h_n) = Delta h_{n+1} - Delta h_n)

推广一下，我们还能得到原序列的(p)阶差分序列
[ Delta^p h_0, Delta^p h_1, dots , Delta^p h_n, dots ]
其中(Delta^p h_n = Delta(Delta^{p-1} h_n) = Delta^{p-1} h_{n+1} - Delta^{p-1} h_n)。

特别的(p = 0)时，(Delta^0 h_n = h_n)。

其中，(p)阶差分在第(p)行上（从第(0)行开始）。

考虑一个序列(h)，通项是一个关于(n)的(p)次多项式，即(h_n = a_{p}n^{p} + a_{p-1}n^{p-1} + cdots + a_0n^0)。

那么这个序列的一阶差分$Delta h $通项就变为了
[ egin{aligned}Delta h_n &= h_{n+1} - h_{n} \&= (a_p (n+1)^p + a_{p-1} (n+1)^{p-1} + cdots + a_0 (n+1)^0) - (a_pn^p + a_{p-1}n^{p-1} + cdots + a_0 n^0)end{aligned} ]
观察$ Delta h_n(的)p(次项：)a_p(n+1)^p - a_p n^p(，把)(n+1)^p$用二项式定理展开后
[ egin{aligned}a_p(n+1)^p - a_p n^p &= a_p(n^p + inom{p}{1} n^{p-1} + inom{p}{2}n^{p-2} + cdots + 1) - a_pn^pend{aligned} ]
发现了什么？(n^p)和(n^p)刚好抵消！即
[ a_p(n+1)^p - a_p n^p = inom{p}{1} a_pn^{p-1} + inom{p}{2}a_pn^{p-2} + cdots + 1 ]
换句话说，每做一次差分，原序列的通项多项式的次数至少会减少(1)。

---

好了来看题吧。

题目让我们求(sumlimits_{i = 1} ^ n i^k)，我们不妨设(h_n = sumlimits_{i = 1} ^ n i^k)，那把这样的序列(h)，做一次差分得到的(Delta h)长啥样呢？

根据定义(Delta h_n = h_{n+1} - h_{n} = sumlimits_{i = 1}^{n+1} i^k - sumlimits_{i = 1}^{n} i^k = (n+1)^k)。显然他的差分序列的通项是一个(k)次多项式，那么(h_n)就是一个(k+1)次多项式。

我们不设(h)序列的通项是一个(k+1)次多项式(f(x))，满足(f(n) = sumlimits_{i = 1}^{n} i^{k})，现在我们要求(f(n))（这里是题目的(n)，上面懒得改了。。。）

观察到(k)只有(10^6)，所以我们可以算出(f(1), f(2), f(3), cdots, f(k+2))，然后拉格朗日差值求出(f(n))就可以了。

但拉格朗日的复杂度是(O(k^2))的。。。还要加一些优化。

我们有公式
[ f(x) = sumlimits_{i = 1}^{k+2} y_i prodlimits_{j ot= i} frac{x - x_j}{x_i - x_j} ]
（这里的(x_i = i, y_i = f(x_i))）。

考虑化简后面的乘法。

我们先看分母。(prodlimits_{j ot= i} frac{1}{x_i - x_j})这个显然是((-1)^{k+2-i} frac{1}{(i-1)!(k+2-i)!})，（后半段是(-1,-2,...,-(k+2-i))的积，直接把(-1)提到前面就好了）题目说了对(10^9 + 7)取模，预处理阶乘逆元就好了。

分子(prodlimits_{j ot= i} x - x_j)。写开来后就是((x-1)(x-2)(x-3)cdots(x-(i-1)) cdot (x-(k+2))(x-(k+1))cdots(x-(i+1)))。

维护两个数组(pre[i] = (x-1)(x-2)(x-3)cdots(x-i),suf[i]=(x-(k+2))(x-(k+1))cdots(x-i))就好了。

前面(y_i)的话可以在枚举(i)的时候带着算，用快速幂。

总复杂度(O(k log k))（快速幂还要(log)啊。。。）

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int P=1e9+7;
const int N=1e6+10;
void update(int &x,int y){
    x+=y; if (x>=P) x-=P;
}
int fpow(int a,int b){
    int ret=1; for (;b;b>>=1,a=1ll*a*a%P) if (b&1)ret=1ll*ret*a%P;
    return ret;
}
#define normal(x) (((x)%P+P)%P)
int pre[N],suf[N],ifac[N],inv[N];
int main(){
    int n,K; scanf("%d%d",&n,&K);
    if (K==0){printf("%d
",n);return 0;} // 特判一下吧，很稳。
    ifac[0]=1,inv[1]=1,ifac[1]=1; // 阶乘的定义，注意0!=1
    for (int i=2;i<=K+2;i++){
        inv[i]=1ll*inv[P%i]*(P-P/i)%P;
        ifac[i]=1ll*ifac[i-1]*inv[i]%P;
    }
    pre[0]=1; // 特殊处理i==0的时候
    for (int i=1;i<=K+2;i++) pre[i]=1ll*pre[i-1]*normal(n-i)%P;
    suf[K+3]=1; // 同上，特殊对待
    for (int i=K+2;i>=0;i--) suf[i]=1ll*suf[i+1]*normal(n-i)%P;
    int yi=0,ans=0; for (int i=1;i<=K+2;i++){
        update(yi,fpow(i,K));
        ll tmp=1ll*yi*ifac[i-1]%P*ifac[K+2-i]%P*pre[i-1]%P*suf[i+1]%P;
        if ((K+2-i)%2==0) update(ans,tmp); else update(ans,P-tmp);
    }
    printf("%d
",ans);
    return 0;
}