题意:求$\begin{align*}\left(\sum\limits_{i=1}^n\dfrac 1i\right)\%\ p^k\end{align*}$
数学真是太可爱了==
直接推公式
设$\begin{align*}f(n,k)=\left(\sum\limits_{i=1}^n\dfrac 1i\right)\%\ p^k\end{align*}$,$\begin{align*}g(n,k)=\left(\sum\limits_{\substack{1\leq i\leq n\\p\nmid i}}\dfrac 1i\right)\%\ p^k\end{align*}$,则有$\begin{align*}f(n,k)=\left(g(n,k)+\dfrac{f\left(\left\lfloor\frac np\right\rfloor,k+1\right)}p\right)\%\ p^k\end{align*}$,下面我们来算$g$
$\begin{align*}g(n,k)&=\sum\limits_{\substack{i=a+bp\\i\leq n}}\dfrac 1i\\&=\sum\limits_{a=1}^{p-1}\sum\limits_{b=0}^{\left\lfloor\frac{n-a}p\right\rfloor}\dfrac 1{a+bp}\\&=\sum\limits_{a=1}^{p-1}\sum\limits_{b=0}^{\left\lfloor\frac{n-a}p\right\rfloor}\dfrac 1a\sum\limits_{i=0}^{k-1}\left(-\dfrac{bp}a\right)^i\\&=\sum\limits_{i=0}^{k-1}(-p)^i\sum\limits_{a=1}^{p-1}\dfrac1{a^{i+1}}\sum\limits_{b=0}^{\left\lfloor\frac{n-a}p\right\rfloor}b^i\\&=\sum\limits_{i=0}^{k-1}(-p)^i\sum\limits_{a=1}^{p-1}\dfrac1{a^{i+1}}S_i\left(\left\lfloor\dfrac{n-a}p\right\rfloor\right)\end{align*}$
最后这里$\begin{align*}S_k(n)=\sum\limits_{i=0}^ni^k\end{align*}$是自然数幂求和,直接算复杂度太高,我们可以推导一番(约定以下出现的第一类斯特林数$\left[\begin{matrix}n\\k\end{matrix}\right]$都是带正负号的)
注意到$\begin{align*}x^{\underline n}=\sum\limits_{j=1}^n\left[\begin{matrix}n\\j\end{matrix}\right]x^j\end{align*}$的$j=n$那一项等于$x^n$,我们可以把它拆出来,把$x^{\underline n}$写成$\begin{align*}n!\binom xn\end{align*}$后移项,我们得到这样的形式:$\begin{align*}x^n=n!\binom xn-\sum\limits_{j=1}^{n-1}\left[\begin{matrix}n\\j\end{matrix}\right]x^j\end{align*}$
引理$\begin{align*}\sum\limits_{l=k}^n\binom lk=\binom{n+1}{k+1}\end{align*}$可用归纳法证
假设当$n=m$时成立,那么$\begin{align*}\sum\limits_{l=k}^{m+1}\binom lk=\binom{m+1}{k+1}+\binom{m+1}k=\binom{m+2}{k+1}\end{align*}$,由归纳法,引理得证
$\begin{align*}S_k(n)&=\sum\limits_{i=0}^n\left(k!\binom ik-\sum\limits_{j=1}^{k-1}\left[\begin{matrix}k\\j\end{matrix}\right]i^j\right)\\&=k!\left(\sum\limits_{i=k}^n\binom ik\right)-\sum\limits_{i=0}^n\sum\limits_{j=1}^{k-1}\left[\begin{matrix}k\\j\end{matrix}\right]i^j\\&=k!\binom{n+1}{k+1}-\sum\limits_{j=1}^{k-1}\left[\begin{matrix}k\\j\end{matrix}\right]\sum\limits_{i=0}^ni^j\\&=\dfrac{(n+1)^{\underline{k+1}}}{k+1}-\sum\limits_{j=1}^{k-1}\left[\begin{matrix}k\\j\end{matrix}\right]S_j(n)\end{align*}$
所以我们$O(k^2)$预处理第一类斯特林数并$O(k^2)$递推求得自然数幂求和
最后看回原式:$\begin{align*}g(n,k)=\sum\limits_{i=0}^{k-1}(-p)^i\sum\limits_{a=1}^{p-1}\dfrac1{a^{i+1}}S_i\left(\left\lfloor\dfrac{n-a}p\right\rfloor\right)\end{align*}$,其中的$\dfrac1{a^{i+1}}$是关于$a$的积性函数,可以线性筛预处理求得,而因为$1\leq a\leq p-1$,所以$\left\lfloor\dfrac{n-a}p\right\rfloor$最多只有两种不同的取值,预处理两次就好了,每次计算$g$的时间复杂度是$O(kp)$,总时间复杂度是$O\left(kp\log_pn\right)$
#include<stdio.h>
typedef long long ll;
ll pow(ll a,ll b){
ll s=1;
while(b){
if(b&1)s*=a;
a*=a;
b>>=1;
}
return s;
}
ll mul(ll a,ll b,ll p){
ll tmp=(ll)((double)a*b/p+1e-6)*p;
return a*b-tmp;
}
ll pow(ll a,ll b,ll p){
ll s=1;
while(b){
if(b&1)s=mul(s,a,p);
a=mul(a,a,p);
b>>=1;
}
return s;
}
ll p,st[100][100],s1[100],s2[100],pw[100010],pr[100010];
bool np[100010];
void stir(ll n,ll k){
ll i,j,mo;
mo=pow(p,k);
st[0][0]=1;
for(i=1;i<=n;i++){
for(j=1;j<=i;j++)st[i][j]=(mul(st[i-1][j],i-1,mo)+st[i-1][j-1])%mo;
}
for(i=1;i<=n;i++){
for(j=1;j<=i;j++){
if((i+j)&1)st[i][j]=-st[i][j];
}
}
}
void gets(ll*s,ll n,ll K){
s[0]=n+1;
ll k,j,mo;
mo=pow(p,K);
for(k=1;k<=K;k++){
s[k]=1;
for(j=n+1;j>n-k;j--)s[k]=mul(s[k],(j/(k+1)*(k+1)==j)?j/(k+1):j,mo);
for(j=1;j<k;j++)s[k]=(s[k]-mul(st[k][j],s[j],mo))%mo;
}
}
void sieve(ll t,ll n,ll mo){
ll M,i,j;
M=0;
pw[1]=1;
for(i=2;i<p;i++){
if(!np[i]){
M++;
pr[M]=i;
pw[i]=pow(pow(i,n,mo),t,mo);
}
for(j=1;j<=M&&i*pr[j]<p;j++){
np[i*pr[j]]=1;
pw[i*pr[j]]=mul(pw[i],pw[pr[j]],mo);
if(i%pr[j]==0)break;
}
}
}
ll g(ll n,ll k){
stir(k,k);
gets(s1,(n-1)/p,k);
if((n-1)/p>0)gets(s2,(n-1)/p-1,k);
ll r,i,a,s,tmp,b,t,mo;
mo=pow(p,k);
r=(n-1)/p;
s=0;
b=1;
t=pow(p,k)-pow(p,k-1)-1;
for(i=0;i<k;i++){
sieve(t,i+1,mo);
tmp=0;
for(a=1;a<p&&n-a>=0;a++){
tmp=(tmp+mul(pw[a],(((n-a)/p==r)?s1:s2)[i],mo))%mo;
}
s=(s+mul(b,tmp,mo))%mo;
b*=-p;
}
return s;
}
ll f(ll n,ll k){
if(n==0)return 0;
return(g(n,k)+f(n/p,k+1)/p)%pow(p,k);
}
int main(){
ll k,n,mo;
scanf("%lld%lld%lld",&p,&k,&n);
mo=pow(p,k);
printf("%lld\n",(f(n,k)+mo)%mo);
}