解决自然数幂和的各种方法

前言

自然数幂和，

暴力

咳咳这个就不用说了。当然不要太暴力了，像求幂可以用快速幂。

高斯消元

我们仔细观察发现，对于求自然数k次方和的公式应该是k+1次的。那么我们可以想到，列个方程组让高斯消元解决，然后得到各项系数后直接代入即可。
复杂度$O(k^3)$

矩阵乘法

发现方法有点问题，先留坑

倍增

我们设$f(n,k)=\sum_{i=1}^ni^k$
那么现在我们要计算f(n,k)怎么办呢？
我们采用分治思想。
如果n是奇数那么
$f(n,k)=f(n-1,k)+n^k$
否则，我们可以先求出n/2的f，然后对于n/2+1~n中的每个数都可以表示为n/2+i
那么
$f(n,k)=f(n/2,k)+\sum_{i=1}^{n/2}\sum_{j=0}^kC_k^j*i^j*(n/2)^{k-j}$
调换一下
$f(n,k)=f(n/2,k)+\sum_{j=0}^kC_k^j*f(n/2,j)*(n/2)^{k-j}$
于是对于每一个n的f我们需要$O(k^2)$，一共有log n层，算上快速幂的复杂度
总共是$O(k^2 log n log k)$

插值法

分拉格朗日插值法和牛顿插值法，我不会……

伯努利数

同不会

第一类斯特林数

该方法来自GEOTCBRL小学生数学题题解，%%%。
我们设

$$Sk(n)=\sum_{i=0}^ni^k$$
那么这个东东要怎么搞呢？
先引入第一类斯特林数的定义：
我们记 $f(n,k)=n*(n-1)*(n-2)*...*(n-k+1)$
$$f(x,n)=\sum_{k=0}^ns_s(n,k)x^k$$
那么我们很容易推出其递推式
$s_s(n,m)=s_s(n-1,m-1)-(n-1)*s_s(n-1,m)$
提一个有趣的一点：如果令 $s_u(n,m)=|s_s(n,m)|$
那么 $s_u(n,m)=s_u(n-1,m-1)+(n-1)*s_u(n-1,m)$
而且有一个性质 $s_u(n,m)=(-1)^{n+m}s_s(n,m)$
而且 $s_u$还有另外的含义， $s_u(n,m)$表示将 n 个不同元素构成m个圆排列的数目。
圆排列是什么？可以理解为两个排列经过旋转可以相同那么视为等价。
由于标程使用了 $s_u$，接下来全程都介绍如何用 $s_u$求自然数幂和。
基本性质：
$s_u(n,n)=1$
$s_u(n,0)=0$
有一个结论：
$k!*C_j^k=f(j,k)$
怎么证明呢？
左边是 $k!*\frac{j!}{k!(j-k)!}=\frac{j!}{(j-k)!}=(j-k+1)*...j=f(j,k)$
那这个结论就显然了。
上面的东西有啥用？等等就会用到。
我们知道 $s_s(k,k)=1$。
那么 $j^k=s_s(k,k)*j^k$
于是呢
$$j^k=k!*C_j^k-f(j,k)+s_s(k,k)*j^k$$
$$j^k=k!*C_j^k-(f(j,k)-s_s(k,k)*j^k)$$
观察括号内的东西
$$f(j,k)-s_s(k,k)*j^k=\sum_{i=0}^ks_s(k,i)*j^i-s_s(k,k)*j^k$$
咦！
$$f(j,k)-s_s(k,k)*j^k=\sum_{i=0}^{k-1}s_s(k,i)*j^i$$
那么
$$j^k=k!*C_j^k-\sum_{i=0}^{k-1}s_s(k,i)*j^i$$
于是
$$Sk(n)=\sum_{j=0}^n(k!*C_j^k-\sum_{i=0}^{k-1}(-1)^{i+k}s_u(k,i)j^i)$$
$$Sk(n)=k!\sum_{j=0}^nC_j^k-\sum_{i=0}^{k-1}(-1)^{i+k}s_u(k,i)\sum_{j=0}^nj^i$$
运用上面的结论
$$Sk(n)=k!C_{n+1}^{k+1}-\sum_{i=0}^{k-1}(-1)^{i+k}s_u(k,i)\sum_{j=0}^nj^i$$
$$Sk(n)=\frac{f(n+1,k+1)}{k+1}-\sum_{i=0}^{k-1}(-1)^{i+k}s_u(k,i)\sum_{j=0}^nj^i$$
我们预处理 $s_u$需要 $O(k^2)$，然后可以用 $O(k^2)$计算 $Sk$。对于减号前面的东西我们直接 $O(k)$做就好了。
如果要求模，该怎么办？
相对于某些方法需要对模数分解质因数，对每一个p^k做了后用中国剩余定理，该方法可以直接对减号前部分暴力扫一遍相乘并模。分母是k+1,分子是k+1个连续自然数相乘，那么一定有一个自然数是k+1的倍数，因此消除分母的影响，让该方法可以轻松解决任意模数的情况。而某些方法要求模数是质数，还有些方法要求模数拆分后最大的质数幂不能太大。而该方法则比较优秀。代码实现方面也比较简单，时间复杂度为 $O(k^2)$，因此个人认为十分好用。
下面是该题做WYF的盒子的代码

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
ll su[2000+10][2000+10],f[2000+10];
ll i,j,k,l,t,n,m,p;
ll qsc(ll x,ll y){
    ll a1=x/1000000,a2=x%1000000,b1=y/1000000,b2=y%1000000;
    ll t=a1*b1%p*1000000%p*1000000%p;
    t=(t+a2*b2%p);
    t=(t+a1*b2%p*1000000%p);
    t=(t+a2*b1%p*1000000%p);
    return t%p;
}
ll qsm(ll x,ll y){
    if (!y) return 1;
    ll t=qsm(x,y/2);
    t=qsc(t,t);
    if (y%2) t=qsc(t,x);
    return t;
}
ll get(ll n){
    if (!n) return 0;
    if (n<=k){
        t=0;
        fo(i,1,n){
            l=1;
            fo(j,1,k)
                l=qsc(l,i);
            t=(t+l)%p;
        }
        return t;
    }
    fo(i,0,k) f[i]=0;
    if (n%2==0) f[1]=qsc(n/2,n+1);
    else f[1]=qsc(n,(n+1)/2);
    f[1]=(f[1]+1)%p;
    f[0]=n%p;
    fo(i,2,k){
        t=1;
        fo(j,n+1-i,n+1)
            if (j%(i+1)==0) t=qsc(t,j/(i+1));else t=qsc(t,j);
        f[i]=t;
        fo(j,0,i-1){
            if ((j+i)%2==0) l=1;else l=-1;
            f[i]=((f[i]-qsc(l*su[i][j],f[j]))%p+p);
        }
        f[i]=f[i]%p;
    }
    return ((f[k]-1)%p+p)%p;
}
int main(){
    scanf("%lld%lld%lld%lld",&k,&n,&m,&p);
    if (n-m<=5000){
        t=0;
        fo(i,m,n) t=(t+qsm(i,k));
        printf("%lld\n",t%p);
        return 0;
    }
    su[0][0]=1;
    fo(i,1,k) su[i][0]=0,su[i][i]=1;
    fo(i,1,k)
        fo(j,1,i-1)
            su[i][j]=(su[i-1][j-1]+qsc(i-1,su[i-1][j]));
    printf("%lld\n",((get(n)-get(m-1))%p+p)%p);
}