题意：求(sum_{i=1}^ni^{-1} mod p^k)的值。保证(p)为质数，(n imes p^k≤10^{18})
题解：神题！细节巨多！
设(f(n,k)=sum_{i=1}^{n}i^{-1} mod p^k)，(g(n,k)=sum_{i=1,i
eq jp}^{n}i^{-1} mod p^k)
则(f(n,k)=g(n,k)+frac{f(lfloorfrac{n}{p}
floor,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}^{lfloorfrac{n-a}{p}
floor}frac{1}{a+bp})
(=sum_{a=1}^{p-1}sum_{b=0}^{lfloorfrac{n-a}{p}
floor}(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)^isum_{a=1}^{p-1}frac{1}{a^{i+1}}sum_{b=0}^{lfloorfrac{n-a}{p}
floor}b^i)
我们令(S_k(n))表示(sum_{i=0}^{n}i^k)
则(g(n,k)=sum_{i=0}^{k-1}(-p)^isum_{a=1}^{p-1}frac{1}{a^{i+1}}S_i(lfloorfrac{n-a}{p}
floor))
如何求(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 imes 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)的项系数为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)^isum_{a=1}^{p-1}frac{1}{a^{i+1}}S_i(lfloorfrac{n-a}{p}
floor))
对于一些零碎的部分，我是直接暴力计算的。
这一题数域可能很大，我写了(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
",f(n,k));
return 0;
}