题面
分析
设$$F(n,k)equivsum_{i=1}nfrac{1}{i}pmod{pk}$$
(p)的倍数都没有逆元，因此一定是把(p)的倍数的倒数里的(p)提出后剩余部分为(p)的倍数。因为题目保证有解，因此我们按这个思路分开求就好了。
不妨先设(n)为(p)的倍数，剩余暴力即可。
再设出不为(p)的倍数的部分$$G(n,k)equivsum_{i=1}{p-1}sum_{j=0}{frac{n}{p}-1}frac{1}{i+jp}pmod{p^k}$$
则有：
[egin{align*}
F(n,k)&equivsum_{i=1}^nfrac{1}{i}\
&equivsum_{i=1}^frac{n}{p}frac{1}{i p}+G(n,k)\
&equivfrac{1}{p}sum_{i=1}^frac{n}{p}frac{1}{i}+G(n,k)\
&equivfrac{F(frac{n}{p},k)}{p}+G(n,k)\
end{align*}]
题目已保证除(p)必定整除。但为了除到(p^k)范围内，将其模数从(p^k)改为(p^{k+1})：$$F(n,k)equivfrac{F(frac{n}{p},k+1)}{p}+G(n,k)$$
第一部分递归即可。则问题转换为如何快速求(G(n,k))。
根据$$frac{1}{1-x}=sum_{k=0}^infty x^k$$
考虑展开(frac{1}{a+bp})：
[egin{align*}
frac{1}{a+bp}&equivfrac{a^{-1}}{1+a^{-1}bp}\
&equiv a^{-1}sum_{i=0}^infty left(-frac{b p}{a}
ight)^i\
&equiv frac{1}{a}sum_{i=0}^infty (-1)^ifrac{b^i}{a^i}p^ipmod{p^k}\
&equiv frac{1}{a}sum_{i=0}^{k-1} (-1)^ifrac{b^i}{a^i}p^i\
end{align*}]
故
[egin{align*}
G(n,k)&equivsum_{i=1}^{p-1}sum_{j=0}^{frac{n}{p}-1}frac{1}{i}sum_{w=0}^{k-1} (-1)^wfrac{j^w}{i^w}p^w\
&equivsum_{i=1}^{p-1}sum_{w=0}^{k-1}(-1)^w p^wfrac{1}{i^{w+1}}sum_{j=0}^{frac{n}{p}-1}j^w
end{align*}]
实验后我们发现应该约定(0^0=1)。
那么我们预处理逆元和(p^w)和自然数幂求和，就可以(O(k p))计算了。
自然数幂求和可以(O(k^2))求出，具体可以看我的博客自然数幂求和(注意边界)
时间复杂度(O(k plog n))
代码
#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
ll n,p,k,m=1,S2[101][101],inv[100001],Snm[101];
inline ll add(ll a,ll b,ll m){return a+b>=m?a+b-m:a+b;}
inline ll mul(ll a,ll b,ll m){return ((a*b-(ll)((double)a/m*b+0.5)*m)%m+m)%m;}
inline ll fp(ll a,ll m){return a&1?m-1:1;}
void init(ll k,ll m){
S2[0][0]=1;
for(int i=1;i<=k;i++)for(int j=1;j<=k;j++)S2[i][j]=add(S2[i-1][j-1],mul(S2[i-1][j],j,m),m);
}
ll S(ll w,ll n,ll m){
if(!w)return n+1;
ll ans=0,facpw=n;
for(int i=1;i<=w;i++){
ans=add(ans,mul(S2[w][i],mul(facpw,inv[i+1],m),m),m);
facpw=mul(facpw,n+m-i,m);
}
return mul(ans,n+1,m);
}
ll g(ll n,ll k,ll m){
ll r=n%p,c=n/p,ans=0;
for(int i=2;i<p;i++)inv[i]=mul(m-m/i,inv[m%i],m);
for(int a=1;a<=r;a++){
ll cnt=0,tpow=1,cc=mul(mul(inv[a],c,m),p,m);
for(int i=0;i<k;i++){
cnt=add(cnt,mul(fp(i,m),tpow,m),m);
tpow=mul(tpow,cc,m);
}
ans=add(ans,mul(cnt,inv[a],m),m);
}
if(c){
init(k-1,m);
for(int w=0;w<k;w++)Snm[w]=S(w,c-1,m);
for(int i=1;i<p;i++){
ll tpow=inv[i],cc=mul(inv[i],p,m);
for(int w=0;w<k;w++){
ans=add(ans,mul(fp(w,m),mul(tpow,Snm[w],m),m),m);
tpow=mul(tpow,cc,m);
}
}
}
return ans;
}
ll f(ll n,ll k,ll m){
if(!n)return 0;
else return add(g(n,k,m),f(n/p,k+1,m*p)/p,m);
}
int main(){
scanf("%lld%lld%lld",&p,&k,&n);
for(int i=1;i<=k;i++)m*=p;
inv[1]=1;
printf("%lld
",f(n,k,m));
}