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

题目和标签都是假的。
题目大意：给定p, k, n$p,k,n$，求：
(∑ni=11i)mod  pk$\left(\sum _{i=1}^n\frac{1}{i}\right)mod{p}^k$
保证有解，pkn≤1018, p≤105$p^{k}n\le {10}^{18},p\le {10}^5$，p$p$是奇素数。
例如，当p=3, k=1, n=6$p=3,k=1,n=6$的时候，尽管13$\frac{1}{3}$和16$\frac{1}{6}$都没有意义，但是他们的和是12$\frac{1}{2}$，就可以计算。
题解：来一波充斥着神奇的推导。

f(n,k)=(∑i=1n1i)mod  pk$$f(n,k)=\left(\sum _{i=1}^n\frac{1}{i}\right)mod{p}^k$$

g(n,k)=(∑i=1,i mod p≠0n1i)mod  pk$$g(n,k)=\left(\sum _{i=1,imodp\ne 0}^n\frac{1}{i}\right)mod{p}^k$$

显然有：

f(n,k)=⎛⎝⎜g(n,k)+f(⌊np⌋,k+1)p⎞⎠⎟mod  pk$$f(n,k)=\left(g(n,k)+\frac{f\left(\lfloor \frac{n}{p}\rfloor ,k+1\right)}{p}\right)mod{p}^k$$

考虑g怎么求：

g(n,k)=(∑i=1,i mod p≠0n1i)mod  pk$$g(n,k)=\left(\sum _{i=1,imodp\ne 0}^n\frac{1}{i}\right)mod{p}^k$$

我们假定i=bp+a, a∈[1,p)$i=bp+a,a\in [1,p)$

g(n,k)=⎛⎝⎜⎜∑a=1p−1∑b=0⌊n−ap⌋1a+bp⎞⎠⎟⎟mod  pk$$g(n,k)=\left(\sum _{a=1}^{p-1}\sum _{b=0}^{\lfloor \frac{n-a}{p}\rfloor }\frac{1}{a+bp}\right)mod{p}^k$$

g(n,k)=⎛⎝⎜⎜∑a=1p−11a∑b=0⌊n−ap⌋11+bpa⎞⎠⎟⎟mod  pk$$g(n,k)=\left(\sum _{a=1}^{p-1}\frac{1}{a}\sum _{b=0}^{\lfloor \frac{n-a}{p}\rfloor }\frac{1}{1+\frac{bp}{a}}\right)mod{p}^k$$

我们注意到：

11−z=∑izi$$\frac{1}{1-z}=\sum _i{z}^i$$

那么令

z=−bpa$$z=-\frac{bp}{a}$$

并代入上式，得到：

g(n,k)=⎛⎝⎜⎜∑a=1p−11a∑b=0⌊n−ap⌋∑i=0k−1(−bpa)i⎞⎠⎟⎟mod  pk$$g(n,k)=\left(\sum _{a=1}^{p-1}\frac{1}{a}\sum _{b=0}^{\lfloor \frac{n-a}{p}\rfloor }\sum _{i=0}^{k-1}{\left(-\frac{bp}{a}\right)}^i\right)mod{p}^k$$

这里内层的i$i$之所以枚举到k−1$k-1$是由于当k≤i$k\le i$时在对pk$p^k$取模结果是0。
继续化简，将i$i$放到外层枚举：

g(n,k)=⎛⎝⎜⎜∑i=0k−1(−p)i∑a=1p−11ai+1∑b=0⌊n−ap⌋bi⎞⎠⎟⎟mod  pk$$g(n,k)=\left(\sum _{i=0}^{k-1}(-p{)}^i\sum _{a=1}^{p-1}\frac{1}{a^{i+1}}\sum _{b=0}^{\lfloor \frac{n-a}{p}\rfloor }{b}^i\right)mod{p}^k$$

g(n,k)=(∑i=0k−1(−p)i∑a=1p−11ai+1Si(⌊n−ap⌋))mod  pk$$g(n,k)=\left(\sum _{i=0}^{k-1}(-p{)}^i\sum _{a=1}^{p-1}\frac{1}{a^{i+1}}{S}_i\left(\lfloor \frac{n-a}{p}\rfloor \right)\right)mod{p}^k$$

其中Sk(n)$S_k(n)$表示从0到n的k次幂和，我们推导这个。注意到因为a∈(0,p)$a\in (0,p)$所以⌊n−ap⌋$\lfloor \frac{n-a}{p}\rfloor$最小值和最大值之差不会超过1。
我们知道xk––$x^{\underset{\_}{k}}$可以用带符号第一类斯特林数展开，即：

xk––=∑i=1ks(k,i)xi$$x^{\underset{\_}{k}}=\sum _{i=1}^{k}s(k,i)x^i$$

那么我们把i=k$i=k$的那一项单独拿出来，就会知道：

xk=xk––−∑i=1k−1s(k,i)xi$$x^k=x^{\underset{\_}{k}}-\sum _{i=1}^{k-1}s(k,i)x^i$$

因此：

Sk(n)=∑x=0nxk=∑x=0n(xk––−∑i=1k−1s(k,i)xi)$$S_k(n)=\sum _{x=0}^n{x}^k=\sum _{x=0}^n\left(x^{\underset{\_}{k}}-\sum _{i=1}^{k-1}s(k,i)x^i\right)$$

我们注意到：

∑x=0nxk––=∑x=0n(xk)k!=1k!∑x=kn(xk)=1k!(n+1k+1)=(n+1)k+1–––––k+1$$\sum _{x=0}^n{x}^{\underset{\_}{k}}=\sum _{x=0}^n\frac{(\genfrac{}{}{0ex}{}{x}{k})}{k!}=\frac{1}{k!}\sum _{x=k}^n(\genfrac{}{}{0ex}{}{x}{k})=\frac{1}{k!}(\genfrac{}{}{0ex}{}{n+1}{k+1})=\frac{(n+1{)}^{\underset{\_}{k+1}}}{k+1}$$

最后一步的转化可以显然的归纳证明是对的。
而：

∑x=0n∑i=1k−1s(k,i)xi=∑i=1k−1s(k,i)∑x=0nxi=∑i=1k−1s(k,i)Si(n)$$\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}^n{x}^i=\sum _{i=1}^{k-1}s(k,i)S_i(n)$$

因此：

Sk(n)=(n+1)k+1–––––k+1−∑i=1k−1s(k,i)Si(n)$$S_k(n)=\frac{(n+1{)}^{\underset{\_}{k+1}}}{k+1}-\sum _{i=1}^{k-1}s(k,i)S_i(n)$$

这样Si(n)(i=0...k)$S_i(n)(i=\mathrm{0...}k)$就可以在O(k2)$O(k^2)$时间内求出。
前面的下降幂下面的除数不能直接逆元算，要特判。
乘法取模不能直接乘，可能会爆炸longlong，因此我开了__int128。
在算g函数的时候需要预处理幂和，这样复杂度显然就是O(pklgpn)$O(pkl{g}_{p}n)$的了。

#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));
}