【GDKOI2016】小学生数学题 【自然数幂求和】【斯特林数】

题意：求$\sum_{i=1}^ni^{-1}\ mod\ p^k$的值。保证$p$为质数，$n\times p^k≤10^{18}$
题解：神题！细节巨多！
设$f(n,k)=\sum_{i=1}^{n}i^{-1}\ mod\ p^k$，$g(n,k)=\sum_{i=1,i\neq jp}^{n}i^{-1}\ mod\ p^k$
则$f(n,k)=g(n,k)+\frac{f(\lfloor\frac{n}{p}\rfloor,k+1)}{p}$
为什么后面那里是k+1呢？
因为若$a≡b\ mod\ c$，则$ak≡bk\ mod\ ck$
我们把那一部分先整体乘$p$，最后再除回来就好。
这样我们就可以递归计算了。
如何求g(n,k)？
$g(n,k)$
$=\sum_{i=a+bp<=n}\frac{1}{i}$
$=\sum_{a=1}^{p-1}\sum_{b=0}^{\lfloor\frac{n-a}{p}\rfloor}\frac{1}{a+bp}$
$=\sum_{a=1}^{p-1}\sum_{b=0}^{\lfloor\frac{n-a}{p}\rfloor}(a+bp)^{-1}$
这个东西用广义二项式定理展开一下
$(a+bp)^{-1}=\sum_{i=0}^{\infty}C_{-1}^{i}a^{-1-i}b^ip^i=\sum_{i=0}^{\infty}(-1)^ia^{-1-i}b^ip^i$
又因为我们模了一个$p^k$
所以
$(a+bp)^{-1}=\sum_{i=0}^{k-1}(-1)^ia^{-1-i}b^ip^i$
所以
$g(n,k)$
$=\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$
我们令$S_k(n)$表示$\sum_{i=0}^{n}i^k$
则$g(n,k)=\sum_{i=0}^{k-1}(-p)^i\sum_{a=1}^{p-1}\frac{1}{a^{i+1}}S_i(\lfloor\frac{n-a}{p}\rfloor)$
如何求$S_k(n)?$
我们先证一个东西：
$x^{\overline n}=\Pi_{i=x}^{x+n-1}i=\sum_{k}s(n,k)x^k$，s是第一类斯特林数。
可以转化为证$[x^p]x^{\overline n}=[x^p]\sum_{k}s(n,k)x^{k}=s(n,p)$
显然在$n=1$的情况下是成立的。
$[x^k]x^{\overline n}$
$=[x^k](x+n-1)x^{\overline{n-1}}$
$=[x^k](x+n-1)x^{\overline{n-1}}$
$=[x^k]x\times x^{\overline{n-1}}+[x^k](n-1)x^{\overline{n-1}}$
$=[x^{k-1}]x^{\overline{n-1}}+[x^k](n-1)x^{\overline{n-1}}$
$=s(n-1,k-1)+(n-1)s(n-1,k)$
$=s(n,k)$
证明完毕。
我们再证一个东西：

$$x^n=n!C_{x}^{n}-\sum_{j=1}^{n-1}s(n,j)x^j$$
其实挺显然的，上面那个证明的式子 $x^n$的项系数为1，所以一减就可以得到这个式子。
我们还要证一个东西：
$$\sum_{i=k}^nC_{i}^{k}=C_{n+1}^{k+1}$$
怎么证？
$$\sum_{i=k}^nC_{i}^{k}=\sum_{i=k}^{n-1}C_{i}^{k}+C_{n}^{k}\\=C_{n}^{k}+C_{n}^{k+1}\\=C_{n+1}^{k+1}$$
证明完毕。
终于要继续自然数幂求和啦！有了之前推的一堆式子，我们就可以化简式子了。
$S_{k}(n)$
$=\sum_{i=0}^{n}i^k$
$=\sum_{i=0}^{n}k!C_{i}^{k}-\sum_{j=1}^{k-1}s(k,j)i^j$
$=k!(\sum_{i=0}^{n}C_{i}^{k})-\sum_{i=0}^{n}\sum_{j=1}^{k-1}s(k,j)i^j$
$=k!C_{n+1}^{k+1}\sum_{j=1}^{k-1}s(k,j)\sum_{i=0}^{n}i^j$
对于固定的 $n$，这个可以$O(k^2)$预处理出来。然后代会这个式子里计算。
$g(n,k)=\sum_{i=0}^{k-1}(-p)^i\sum_{a=1}^{p-1}\frac{1}{a^{i+1}}S_i(\lfloor\frac{n-a}{p}\rfloor)$
对于一些零碎的部分，我是直接暴力计算的。
这一题数域可能很大，我写了 $O(1)$快速乘法。
这题好难写啊！QAQ
代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define int long long
using namespace std;
int p,k,n,x,y,r[100005],s[105][105];
int mul(int a,int b,int mod){
    int c=a*(double)b/mod+0.5;
    int res=a*b-c*mod;
    if(res<0){
        res+=mod;
    }
    return res;
}
int fastpow(int a,int x,int mod){
    a%=mod;
    int res=1;
    while(x){
        if(x&1){
            res=mul(res,a,mod);
        }
        x>>=1;
        a=mul(a,a,mod);
    }
    return res;
}
int exgcd(int a,int b,int &x,int &y){
    if(!b){
        x=1;
        y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
int getinv(int a,int b){
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
void init(int n,int k){
    const int mod=pow(p,k);
    r[0]=n+1;
    s[0][0]=1;
    for(int i=1;i<=k;i++){
        for(int j=1;j<=k;j++){
            s[i][j]=(s[i-1][j-1]+mul(s[i-1][j],i-1,mod))%mod;
        }
    }
    for(int i=1;i<=k;i++){
        for(int j=1;j<=i;j++){
            if((i+j)&1){
                s[i][j]=mod-s[i][j];
            }
        }
    }
    for(int i=1;i<=k;i++){
        r[i]=1;
        bool flag=false;
        for(int j=n+1;j>=n-i+1;j--){
            if(!flag&&j%(i+1)==0){
                flag=true;
                r[i]=mul(r[i],j/(i+1),mod);
            }else{
                r[i]=mul(r[i],j,mod);
            }
        }
        if(!flag){
            r[i]=mul(r[i],getinv(i+1,mod),mod);
        }
        for(int j=1;j<i;j++){
            r[i]=(r[i]-mul(s[i][j],r[j],mod)+mod)%mod;
        }
    }
    return;
}
int g(int n,int k){
    const int mod=pow(p,k);
    int rem=n%p,res=0;
    for(int i=n-rem+1;i<=n;i++){
        res=(res+getinv(i,mod))%mod;
    }
    n-=rem;
    if(!n){
        return res;
    }
    init((n-1)/p,k);
    for(int i=0,base=1;i<k;i++,base=(base*(-p)%mod+mod)%mod){
        for(int j=1;j<p;j++){
            res=(res+mul(mul(base,getinv(fastpow(j,i+1,mod),mod),mod),r[i],mod))%mod;
        }
    }
    return res;
}
int f(int n,int k){
    if(!n){
        return 0;
    }
    return (g(n,k)+f(n/p,k+1)/p)%(int)(pow(p,k));
}
signed main(){
    scanf("%lld%lld%lld",&p,&k,&n);
    printf("%lld\n",f(n,k));
    return 0;
}