https://gmoj.net/senior/#main/show/4367
题解：
退役前来了结这个OI入（弃）坑题吧。
(sum_{i=1}^n frac{1}{i} (~mod~p^k))
类似于扩展Lucas，把(p)的倍数和不是(p)的倍数分开。
(= frac{1}{p} sum_{i=1}^{n/p} frac{1}{i} + sum_{i=n/p*p+1}^n frac{1}{i} + sum_{i=0}^{n/p-1}sum_{j=1}^{p-1} frac{1}{ip+j})
这题保证最后分数有逆元，所以第一部分可直接递归((n/p,p^{k+1}))然后除以(p)即可。
第二部分(<p)，所以暴力逆元即可。
现在看第三部分：
(sum_{i=0}^{n/p-1}sum_{j=1}^{p-1} frac{1}{ip+j})
(=sum_{i=0}^{n/p-1}sum_{j=1}^{p-1} j^{-1}frac{1}{ipj^{-1}+1})
(i=0)时，后面(=sum_{j=1}^{p-1} j^-1)，好算。
(sum_{i=1}^{n/p-1}sum_{j=1}^{p-1} j^{-1}frac{1}{ipj^{-1}+1})
生成函数中有(frac{1}{1-x}=sum_{qge 0} x^q)，所以：
(=sum_{i=1}^{n/p-1}sum_{j=1}^{p-1} j^{-1}frac{1}{1-(-ipj^{-1})})
(=sum_{i=1}^{n/p-1}sum_{j=1}^{p-1} j^{-1} sum_{qge 0} (-ipj^{-1})^q)
因为里面有(p)，所以当(q ge k)时，一定是(0)，缩小(q)的循环范围，
(=sum_{i=1}^{n/p-1}sum_{j=1}^{p-1} j^{-1} sum_{q=0}^{k-1} (-ipj^{-1})^q)
(sum_{j=1}^{p-1} j^{-1} sum_{q=0}^{k-1} (-pj^{-1})^q sum_{i=1}^{n/p-1}i^q)
自然数幂和用第二类斯特林数搞搞就好了。
求逆元用exgcd，需要黑科技乘法。
分析时间复杂度：
-
有一个递归，最多递归(log_p n)层。
-
每一层第二部分是(O(kp log p))，第三部分是(O(k^3+pk))。
很明显(O(能过）)。
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i < _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
ll mul(ll x, ll y, ll mo) {
x = (x % mo + mo) % mo, y = (y % mo + mo) % mo;
ll z = (long double) x * y / mo;
z = x * y - z * mo;
if(z < 0) z += mo; else if(z >= mo) z -= mo;
return z;
}
const int N = 505;
ll s[N][N];
ll solve(ll n, ll m, ll mo) {
if(m == 0) return n % mo;
if(n == 0) return 0;
s[0][0] = 1;
fo(i, 1, m) fo(j, 1, i) s[i][j] = (s[i - 1][j - 1] + s[i - 1][j] * j) % mo;
ll ans = 0;
fo(j, 1, m) {
ll xs = s[m][j];
ll n0 = (n + 1) - (n + 1) % (j + 1);
for(ll x = n - j + 1; x < n0; x ++) xs = mul(xs, x, mo);
for(ll x = n0 + 1; x <= n + 1; x ++) xs = mul(xs, x, mo);
xs = mul(xs, n0 / (j + 1), mo);
ans = (ans + xs) % mo;
}
return ans;
}
void exgcd(ll a, ll b, ll &x, ll &y) {
if(!b) { x = a, y = 0; return;}
exgcd(b, a % b, y, x); y -= (a / b) * x;
}
ll ginv(ll a, ll p) {
ll x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
}
ll ksm(ll x, ll y, ll mo) {
ll s = 1;
for(; y; y /= 2, x = mul(x, x, mo))
if(y & 1) s = mul(s, x, mo);
return s;
}
ll dg(ll p, ll k, ll n) {
if(n == 0) return 0;
ll mo = 1; fo(i, 1, k) mo = mo * p;
ll ans = dg(p, k + 1, n / p) / p;
for(ll x = n / p * p + 1; x <= n; x ++) ans = (ans + ginv(x, mo)) % mo;
static ll t[N];
if(n / p > 0) {
// fo(i, 0, n / p - 1) fo(j, 1, p - 1)
// ans = (ans + ginv(i * p + j, mo)) % mo;
fo(q, 0, k - 1) t[q] = solve(n / p - 1, q, mo);
fo(j, 1, p - 1) {
ll invj = ginv(j, mo);
ans = (ans + invj) % mo;
fo(q, 0, k - 1) {
ll xs = mul(invj, ksm(-mul(invj, p, mo), q, mo), mo);
ans = (ans + mul(xs, t[q], mo)) % mo;
}
}
}
ans = (ans % mo + mo) % mo;
return ans;
}
ll p, k, n;
int main() {
scanf("%lld %lld %lld", &p, &k, &n);
ll ans = dg(p, k, n);
pp("%lld
", ans);
}