JZOJ4367. 【GDKOI2016】小学生数学题（数论推理+自然数幂和）

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，需要黑科技乘法。

分析时间复杂度：

1. 有一个递归，最多递归(log_p n)层。

2. 每一层第二部分是(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);
}