高效求解自然数k次幂和的方法

前言

寒假集训比赛时遇到了一道《WYF的盒子》，题目大意就是求$\sum_{i=m}^n i^k$对$p$取模的结果。
其中有的数据点满足$n-m\le 5000$但$k\le{10}^{12}$，剩下的数据点满足$n,m\le10^{12}$且$k\le 2000$。对于所有数据模数$p\le10^{12}$。
其中$n$和$m$相差很小的点显然可以使用快速幂，但是剩下的点怎么做？换言之，我们需要找到一个时间复杂度之和$k$有关的高效的求自然数幂和的算法。

---

普通公式

运用二项式定理我们可以推出一般的公式。

$$S_k=\frac{1}{k+1}[(n+1)^{k+1}-(n+1)-\sum_{i=1}^{k-1}C_{k+1}^{k-i}S_i]$$
推算过程链接：http://wenku.baidu.com/link?url=F07Cj1SxmyNr8_EriB3UrJU8tn4f3aByvpYKPvrLCoqV2HWjeXG1y6XWbgFQL8mBcP5Z0is1Sqpe0zxbDRmMSsBjt-0mQRRvvHTMqq57gpK
当然我们也可以不使用二项式定理，只研究递推关系，也可以得到相同的公式。
推算过程链接：http://blog.csdn.net/acdreamers/article/details/38929067
这种公式固然优美，时间复杂度为$\mathrm O(k^2)$，但是不可避免地要做除法，如果除数和模数互质那好办，直接逆元即可。但是这里显然不满足。

---

拉格朗日插值法&牛顿插值法

暑假时HeD大神给我们讲了这两种解法，初中党表示一脸懵逼。
拉格朗日插值法其实不难，拿高中选修数论初步看看就懂了。但是还是避免不了除法。
牛顿插值法只需要差分，不需要做除，时间复杂度$\mathrm O(k^2)$，是很优秀的算法，然而我不会。

---

伯努利数

这个我也不会，但是我有链接：http://blog.csdn.net/acdreamers/article/details/38929067
然而其实它也要做除法。

---

矩阵乘法

听说可以，然而我并不能YY出来。

---

第一类斯特林数

BB了这么多，终于到重点了。
第一类斯特林数是什么鬼呢？它有一个圆环排列的定义，但是我并不知道。
于是我便引用它最初发现时的定义。
记

$$f(n,m)=\prod_{i=n-m+1}^n i$$
即 $f(n,m)=\frac{n!}{(n-m)!}=n(n-1)(n-2)...(n-m+1)$
则第一类斯特林数就表示 $n$的各次幂的系数。
其实第一类斯特林数分两种：有符号$s_s$和无符号 $s_u$。有符号就是正常的第一类斯特林数，无符号就是有符号的绝对值。
即
$$f(n,m)=\sum_{i=0}^ms_s(m,i)n^i$$
可以发现 $s_s(n,m)=(-1)^{n+m}s_u(n,m)$。不过我们下面不需要使用 $s_u$。
至于 $s_s$和 $s_u$的递推公式，就由大家自己推算了。利用第一类斯特林数的定义很容易推出来，实在不会就参考我的代码。
第一类斯特林数的递推复杂度是 $\mathrm O(k^2)$的。
再讲一个显然的东西。
$$k!C_j^k=f(j,k)$$
这个拆开来约下分就知道了。
还有一个（并不）显然的东西。
$$C_{n+1}^{k+1}=\sum_{i=0}^nC_i^k$$
这个使用组合数的杨辉三角递推式拆开就知道了。

下面开始推算

$$\begin{align} f(j,k)-s_s(k,k)j^k&=\sum_{i=0}^ks_s(k,i)j^i-s_s(k,k)j^k\\ &=\sum_{i=0}^{k-1}s_s(k,i)j^i \end{align}$$
显然 $s_s(k,k)=1$，于是我们移项之后可以得到
$$j^k=f(j,k)-\sum_{i=0}^{k-1}s_s(k,i)j^i$$
对 $j$求和，记$S_k(n)$表示 $1$至$n$的 $k$次幂和得
$$\begin{align} S_k(n)&=\sum_{j=0}^{n}(f(j,k)-\sum_{i=0}^{k-1}s_s(k,i)j^i)\\ &=\sum_{j=0}^{n}(k!C_j^k-\sum_{i=0}^{k-1}s_s(k,i)j^i)\\ &=k!\sum_{j=0}^nC_j^k-\sum_{i=0}^{k-1}s_s(k,i)S_i(n)\\ &=k!C_{n+1}^{k+1}-\sum_{i=0}^{k-1}s_s(k,i)S_i(n)\\ &=\frac{f(n+1,k+1)}{k+1}-\sum_{i=0}^{k-1}s_s(k,i)S_i(n) \end{align}$$
这里面的 $f$可以$\mathrm O(k)$扫一遍，里面一定有且仅有一个 $k+1$的倍数，乘的时候处理一下即可。后面部分 $\mathrm O(k)$枚举，所以总的复杂度为 $\mathrm O(k^2)$。
至此问题完美解决。

---

代码实现

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

typedef long long LL;

const int K=2050;

LL ss[K][K],sum[K];
LL k,n,m,p,ans;

const int A=1000000;

LL mult(LL x,LL y)
{
    LL a1=x/A,a2=x%A,b1=y/A,b2=y%A;
    LL ret=a1*b1%p*A%p*A%p;
    ret+=a1*b2%p*A%p;
    ret+=a2*b1%p*A%p;
    ret+=a2*b2%p;
    return ret%p;
}

LL quick_power(LL x,LL y)
{
    LL ret=1;
    while (y)
    {
        if (y&1)
            ret=mult(ret,x);
        x=mult(x,x);
        y>>=1;
    }
    return ret;
}

void pre()
{
    ss[0][0]=1;
    for (int i=1;i<=k;i++)
        ss[i][0]=0,ss[i][i]=1;
    for (int i=1;i<=k;i++)
        for (int j=1;j<i;j++)
            ss[i][j]=ss[i-1][j-1]+mult(-i+1,ss[i-1][j]);
}

LL calc(LL n)
{
    if (n&1)
        sum[1]=mult(n,(n+1)/2);
    else
        sum[1]=mult(n/2,n+1);
    for (int i=2;i<=k;i++)
    {
        sum[i]=1;
        for (LL j=n+1;j>=n-i+1;j--)
            if (j%(i+1))
                sum[i]=mult(sum[i],j);
            else
                sum[i]=mult(sum[i],j/(i+1));
        for (int j=0;j<i;j++)
            sum[i]=(sum[i]-mult(ss[i][j],sum[j])+p)%p;
    }
    return sum[k];
}

int main()
{
    freopen("box.in","r",stdin);
    freopen("box.out","w",stdout);
    scanf("%lld%lld%lld%lld",&k,&n,&m,&p);
    if (n-m<=5000)
        for (LL i=m;i<=n;i++)
            (ans+=quick_power(i,k))%=p;
    else
    {
        pre();
        ans=(calc(n)-calc(m-1)+p)%p;
    }
    printf("%lld\n",ans);
    fclose(stdin);
    fclose(stdout);
    return 0;
}