小学生数学题 - 用第一类斯特林数解决自然数幂求和这个标签真长

题目和标签都是假的。
题目大意：给定$p,\ k,\ n$，求：
$\left(\sum_{i=1}^n\frac1i\right)mod\ \ p^k$
保证有解，$p^kn\le10^{18},\ p\le10^5$，$p$是奇素数。
例如，当$p=3,\ k=1,\ n=6$的时候，尽管$\frac13$和$\frac16$都没有意义，但是他们的和是$\frac12$，就可以计算。
题解：来一波充斥着神奇的推导。

$$f(n,k)=\left(\sum_{i=1}^n \frac1i\right)mod\ \ p^k$$
$$g(n,k)=\left(\sum_{i=1,i\ mod\ p\neq0}^n\frac1i\right)mod\ \ p^k$$
显然有：
$$f(n,k)=\left(g(n,k)+\frac{f\left(\left\lfloor\frac np\right\rfloor,k+1\right)}p\right)mod\ \ p^k$$
考虑g怎么求：
$$g(n,k)=\left(\sum_{i=1,i\ mod\ p\neq0}^n\frac1i\right)mod\ \ p^k$$
我们假定$i=bp+a,\ a\in[1,p)$
$$g(n,k)=\left(\sum_{a=1}^{p-1}\sum_{b=0}^{\left\lfloor\frac{n-a}p\right\rfloor}\frac{1}{a+bp}\right)mod\ \ p^k$$
$$g(n,k)=\left(\sum_{a=1}^{p-1}\frac1a\sum_{b=0}^{\left\lfloor\frac{n-a}p\right\rfloor}\frac{1}{1+\frac{bp}{a}}\right)mod\ \ p^k$$
我们注意到：
$$\frac{1}{1-z}=\sum_{i}z^i$$
那么令
$$z=-\frac{bp}{a}$$
并代入上式，得到：
$$g(n,k)=\left(\sum_{a=1}^{p-1}\frac1a\sum_{b=0}^{\left\lfloor\frac{n-a}p\right\rfloor}\sum_{i=0}^{k-1}\left(-\frac{bp}{a}\right)^i\right)mod\ \ p^k$$
这里内层的$i$之所以枚举到$k-1$是由于当$k\le i$时在对$p^k$取模结果是0。
继续化简，将$i$放到外层枚举：
$$g(n,k)=\left(\sum_{i=0}^{k-1}(-p)^i\sum_{a=1}^{p-1}\frac{1}{a^{i+1}}\sum_{b=0}^{\left\lfloor\frac{n-a}p\right\rfloor}b^i\right)mod\ \ p^k$$
$$g(n,k)=\left(\sum_{i=0}^{k-1}(-p)^i\sum_{a=1}^{p-1}\frac{1}{a^{i+1}}S_i\left(\left\lfloor\frac{n-a}p\right\rfloor\right)\right)mod\ \ p^k$$
其中$S_k(n)$表示从0到n的k次幂和，我们推导这个。注意到因为$a\in(0,p)$所以$\left\lfloor\frac{n-a}p\right\rfloor$最小值和最大值之差不会超过1。
我们知道$x^{\underline{k}}$可以用带符号第一类斯特林数展开，即：
$$x^{\underline{k}}=\sum_{i=1}^ks(k,i)x^i$$
那么我们把$i=k$的那一项单独拿出来，就会知道：
$$x^k=x^{\underline{k}}-\sum_{i=1}^{k-1}s(k,i)x^i$$
因此：
$$S_k(n)=\sum_{x=0}^n x^k=\sum_{x=0}^n\left(x^{\underline{k}}-\sum_{i=1}^{k-1}s(k,i)x^i\right)$$
我们注意到：
$$\sum_{x=0}^nx^{\underline{k}}=\sum_{x=0}^n\frac{\binom{x}{k}}{k!}=\frac{1}{k!}\sum_{x=k}^n\binom{x}{k}=\frac{1}{k!}\binom{n+1}{k+1}=\frac{(n+1)^{\underline{k+1}}}{k+1}$$
最后一步的转化可以显然的归纳证明是对的。
而：
$$\sum_{x=0}^n\sum_{i=1}^{k-1}s(k,i)x^i=\sum_{i=1}^{k-1}s(k,i)\sum_{x=0}^nx^i=\sum_{i=1}^{k-1}s(k,i)S_i(n)$$
因此：
$$S_k(n)=\frac{(n+1)^{\underline{k+1}}}{k+1}-\sum_{i=1}^{k-1}s(k,i)S_i(n)$$
这样$S_i(n)(i=0...k)$就可以在$O(k^2)$时间内求出。
前面的下降幂下面的除数不能直接逆元算，要特判。
乘法取模不能直接乘，可能会爆炸longlong，因此我开了__int128。
在算g函数的时候需要预处理幂和，这样复杂度显然就是$O(pklg_pn)$的了。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<assert.h>
#define lint long long
#define MAXK 200
#define MAXP 100010
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
lint s[MAXK][MAXK],pk[MAXK],ivk[MAXP][MAXK];
inline lint tms(lint x,lint y,lint k) { return (__int128)x*y%pk[k]; }
inline lint fast_pow(lint x,lint k,lint qk)
{
    lint ans=1ll;
    for(;k;k>>=1,x=tms(x,x,qk))
        if(k&1ll) ans=tms(ans,x,qk);
    return ans;
}
inline lint prelude_stirling(lint n,lint k)
{
    s[1][1]=1;
    for(lint i=2;i<=n;i++)
        for(lint j=1;j<=i;j++)
            s[i][j]=(tms(s[i-1][j],i-1,k)+s[i-1][j-1])%pk[k];
    for(lint i=1;i<=n;i++)
        for(lint j=1;j<=i;j++)
            if((i+j)&1) s[i][j]=(pk[k]-s[i][j])%pk[k];
    return 0ll;
}
lint sk1[MAXK],sk2[MAXK];
lint getS(lint k,lint n,lint qk,lint *sk)
{
    memset(sk,0,sizeof(lint)*(k+1)),sk[0]=(n+1)%pk[qk];
    if(n<=0||k==0) return 0;
    if(n&1) sk[1]=tms((n+1)/2,n,qk);
    else sk[1]=tms(n/2,n+1,qk);
    for(lint i=2;i<=k;i++)
    {
        if(n-i+1<=0) sk[i]=0;
        else{
            sk[i]=1;lint nc=1;
            for(lint j=n+1;j>=n-i+1;j--)
                if(j%(i+1)==0&&nc) sk[i]=tms(sk[i],j/(i+1),qk),nc=0;
                else sk[i]=tms(sk[i],j,qk);
        }
        for(lint j=1;j<i;j++)
        {
            sk[i]=(sk[i]-tms(sk[j],s[i][j]%pk[qk],qk))%pk[qk];
            if(sk[i]<0) sk[i]+=pk[qk];
        }
    }
    return 0ll;
}
lint G(lint n,lint k)
{
    lint ans=0ll,r=n%pk[1];
    getS(k-1,n/pk[1]-1,k,sk1),getS(k-1,n/pk[1],k,sk2);
    for(lint i=0;i<k;i++)
    {
        lint pi=pk[i],s=0ll;if(i&1) pi=(pk[k]-pi)%pk[k];
        for(lint a=1;a<=min(pk[1]-1,n);a++)
            if(a<=r) s=(s+tms(ivk[a][i+1]%pk[k],sk2[i],k))%pk[k];
            else s=(s+tms(ivk[a][i+1]%pk[k],sk1[i],k))%pk[k];
        ans=(ans+tms(pi,s,k))%pk[k];
    }
    return ans;
}
lint F(lint n,lint k) { return n?(G(n,k)+F(n/pk[1],k+1)/pk[1])%pk[k]:0; }
int main()
{
    lint p,k,n;scanf("%lld%lld%lld",&p,&k,&n);
    lint qk=k,ns=n;while(ns) ns/=p,k++;
    for(lint i=pk[0]=1;i<=k;i++) pk[i]=pk[i-1]*p;
    prelude_stirling(k,k);
    for(lint i=1;i<p;i++)
    {
        ivk[i][0]=1,ivk[i][1]=fast_pow(i,pk[k]-pk[k-1]-1,k);
        for(lint j=2;j<=k;j++) ivk[i][j]=tms(ivk[i][j-1],ivk[i][1],k);
    }
    return !printf("%lld\n",F(n,qk));
}