自然数方幂和公式及其快速算法

自然数方幂和公式

关于自然数方幂和公式，网上的求解版本有很多种。这里介绍一种不为人知，十分简洁明了的求解方法，该公式并非原创，但是整个证明过程和方法完全原创。它的思想来源于我高中时在一本数学竞赛书中的数列例题(书名忘了...)，正因为一本本这样的书，让我大学选择了数学系，现在依然在学习数学。

若 \(a_n = n(n-1)\) 求其前 \(n\) 项和 \(S_n\)

\[ a_n = n(n-1) ＝ \frac{(n+1)n(n-1) -n(n-1)(n-2)}{3} \]

所以

\[ S_n =\frac{(n+1)n(n-1)}{3} \]

受到上面做法的启发，我们推广到一般形式：

若 \(a_n = A_n ^p\)，则其前 \(n\) 项和 \(S_n ＝ \frac{ A_{n+1} ^{p+1}}{p+1}\)

\[ a_n = A_n ^p ＝ \frac{A_{n+1} ^{p+1} -A_n ^{p+1}}{p+1} \]

所以

\[ S_n = \frac{A_{n+1} ^{p+1}}{p+1} \]

由上面的结论，我们可以直奔主题了

求 \(1^2+2^2+ \dots + n^2\)

由 \[ n^2 = n(n-1) + n \] 知 \[ 1^2+2^2+ \dots + n^2 = \frac{(n+1)n(n-1)}{3} + \frac{(n+1)n}{2} = \frac{n(n+1)(2n+1)}{6} \]

同理我们可以推广上面做法

求\(1^p+2^p+ \dots + n^p\)

设

\[ n^p = a_p A_n ^p + a_{p-1} A_n ^{p-1} + a_{p-2} A_n ^{p-2} + \cdots + a_1 A_n ^1 \] 则

\[ 1^p+2^p+ \dots + n^p = \frac{a_p A_{n+1} ^{p+1}}{p+1} + \frac{a_{p-1} A_{n+1} ^{p}}{p} + \cdots + \frac{a_1 A_{n+1} ^2}{2} \] 因此问题的关键就转化成如何求解数组 \(a_k ,k=1,2,\cdots,p\) 我们发现当 \(n＝k\) 时成立 \[ k^p = a_k A_k ^k + a_{k-1} A_k ^{k-1} + \cdots a_1 A_{k} ^1 \] 令 \(b_k = k! \cdot a_k\) 则

\[ k^p = b_k + b_{k-1} C_k ^{k-1} + \cdots b_1 C_{k} ^1 \] 即

\[ k^p = \sum _{j=1} ^k b_j C _{k} ^j \] 应用二项式反演 [我的博文](../inverseFormulaAndItsMatrixForm）

\[ b_k = \sum _{j=1} ^k (-1)^{k-j} C _{k} ^j j^p \] 因此最终，我们有公式

\[ 1^p+2^p+ \dots + n^p = \sum _{k=1} ^p \; (\; \sum_{j=1} ^ {k} (-1)^{k-j} C_k^j j^p \;) \; C _{n+1} ^{k+1} \]

妙呀，帅的呀，猛的呀，不谈了呀，哈哈哈

这个东西本质叫 离散微积分，也叫下降幂。

生成函数，离散微积分，上升幂下降幂，Stirling 数，牛顿多项式，特别的有用

自然数方幂和快速算法

这次写的目的是因为上次是从数学上完美的解决了这个问题，这次我们要从计算上完美的解决这个问题，当然这归功于我看到的一份 sgtlaugh 的代码。经过解读体会到其中的奥秘，特此记录。一句话，简直不敢相信。 > 如果有人说他能在 \(O(k)\) 时空复杂度求解 \(\sum_{i=1}^n i^k\)，你肯定会说这怎么可能别忽悠我了，那我只能说，因为你没看过这篇博文。

首先预处理，\(O(k)\) 复杂度 求 \(\sum_{i=1}^n i^k\) 其中 \(n \leq k\)

我之前一直以为要用 \(k \log k\) 的复杂度才能解决这个问题，其实我们只需对所有素数 \(p\) 计算 \(p^k\) 即可。对于一般的 \(i\) 我们先预处理 其最小素因子 \(sp[i]\)。计算 \(sp[i]^k \cdot (i/sp[i])^k\) 即可（具体可见最后代码）。由于素数的阶为 \(O(\frac{k}{\log k})\) 因此整个复杂度即为 \(O(k)\)。

再由 Lagrange 插值多项式得出最终答案

因为我们知道 \(\sum_{i=1} ^n i^k\) 一定是一个关于 \(n\) 的次数为 \(k+1\) 的多项式。因此，我们只需计算其在 \(0,\cdots,k+1\) 上的取值，用 Lagrange 插值多项式即可知道答案。

对于一个次数不超过 \(n\) 的多项式 \(f(x)\)，其在不同位置 \(x_0,\cdots,x_n\) 的取值唯一决定了这个多项式:

\[ f(x) = \sum_{i=0} ^n f(x_i) \prod_{j=0,j \neq i} ^n \frac{x-x_j}{x_i - x_j} \]

具体到本问题，我们取 \(x=n, m=k+1, x_i=i\) 那么

\[ f(n) = \sum_{i=0} ^m (-1)^{m-i} f(x_i) {n \choose i} {n-i-1 \choose m-i } \]

例题：Codeforces 622F

#include <bits/stdc++.h>
typedef long long LL;
using namespace std;
const int N = 1e6+6;
const LL mod = 1e9+7;
int sp[N],p[N];
int spf(){ // samllest prime factor
  int cnt=0;p[cnt++]=2;
  for(int i=2;i<N;i+=2) sp[i]=2;
  for(int i=1;i<N;i+=2) sp[i]=i;
  for(int i=3;i<N;i+=2){
    if(sp[i]==i)    p[cnt++] = i;
    for(int j = 1;j < cnt && p[j]<=sp[i] && i * p[j] < N; ++j) {
      sp[i * p[j]] = p[j];
    }
  }
  return cnt;
}
LL pow_mod(LL x,LL n,LL p){
  LL r=1;
  while(n){
    if(n&1) r=r*x%p;
    n>>=1;  x=x*x%p;
  }
  return r;
}
LL inv[N],AP[N],AS[N],f[N];
LL getpowsum(LL n,int k){ // mod > k
  if(k==0)    return n%mod;
  if(p[0]!=2) spf();
  int nk=k+1;
  f[0]=0;f[1]=1;
  for(int i=2;i<=nk;++i){
    if(sp[i]==i) f[i]=pow_mod(i,k,mod);
    else         f[i]=f[sp[i]]*f[i/sp[i]]%mod;
  }
  for(int i=1;i<=nk;++i){
    f[i]+=f[i-1];
    if(f[i]>=mod) f[i]-=mod;
  }
  if(n<=nk) return f[n];
  LL tmp = 1;
  for(int i=2;i<=nk;++i)    tmp=tmp*i%mod;
  inv[nk] = pow_mod(tmp,mod-2,mod);
  for(int i=nk-1;i>=0;--i) inv[i]=inv[i+1]*(i+1)%mod;
  AP[0]=AS[nk]=1;
  for(int i=1;i<=nk;++i)   AP[i]=AP[i-1]*(n+1-i)%mod;
  for(int i=nk-1;i>=0;--i) AS[i]=AS[i+1]*(n-i-1)%mod;
  LL res = 0;
  for(int i=1;i<=nk;++i){ // because f(i)=0
    LL x = f[i]*AP[i]%mod*AS[i]%mod*inv[i]%mod*inv[nk-i]%mod;
    if((nk-i)&1) res-=x; // be careful
    else        res+=x;
  }
  return (res%mod+mod)%mod;
}
int main(){
  LL n;
  int k;
  while(cin>>n>>k){
    cout<<getpowsum(n,k)<<endl;
  }
  return 0;
}

实际上我们可以不求 \(\mod p\) 后的答案，利用大数类得到标准答案，但是这时因为数字实在太大，每次乘法的用时过大，因此仅适合 \(k<n\) 的情况

cpp_int f[N];
cpp_int getpowsum(LL n,int k){
  if(k==0)    return cpp_int(n);
  if(p[0]!=2) spf();
  int nk=2*k+1;
  f[0]=0;f[1]=1;
  for(int i=2;i<=nk+1;++i){
    if(sp[i]==i) f[i]=pow(cpp_int(i),k);
    else         f[i]=f[sp[i]]*f[i/sp[i]];
  }
  for(int i=1;i<=nk;++i) f[i]+=f[i-1];
  if(n<=nk)   return f[n];
  cpp_int res = 0,tl=1,tr=1;
  for(int i=nk-1;i>=0;--i) tr=tr*(n-i-1)/(nk-i);
  for(int i=0;i<=nk;++i){
    if((nk-i)&1) res -= f[i]*tl*tr;
    else         res += f[i]*tl*tr;
    tl = tl*(n-i)/(i+1);
    tr = tr*(nk-i)/(n-i-1);
  }
  return res;
}

其实如果我们知道最终的上界，求出多个 \(\mod p\) 后的答案，再用中国剩余定理貌似很不错。

该方法可以推广成求 \(\sum_{i=1}^n f(i)^k\)，其中 \(f(x)\) 是多项式。具体分析即可

这种情况一般很难再做到 \(O(k)\) 时间复杂度，而变成了 \(O(k \log k) \deg f\) 复杂度。

如果求所有 \(f(k) = \sum_{i = 1}^n i^k, k = 0, \cdots, m\)，可以在 \(O(m \log m)\) 复杂度求出

\(\frac{f(k)}{k!}\) 是 \(\sum_{i = 1}^n e^{ix} = e^x \frac{e^{nx} - 1}{e^x - 1}\) 的 \(x^k\) 的系数。注意到 \(e^x - 1\) 的常数项系数为 0，所以不可逆，要分子分母同时除以 \(x\)，所以求它们的时候要多预算一位。

这里用幂生成函数而非常规生成函数是因为 \(\sum_{i = 1}^n \frac{1}{1 - i x}\) 不好计算，维护分子分母分治可以在 \(O(n \log^2 n)\) 计算出结果。