CSP-S 模拟 轻飘飘的时间 （莫比乌斯反演）（杜教筛）（Lucas）

题意：给定 $n,m,d$，求
$$\sum _{i=1}^n\sum _{j=1}^m\left(\genfrac{}{}{0ex}{}{gcd(i,j)}{B}\right)$$
直接反演
$$\sum _{d=1}^n\left(\genfrac{}{}{0ex}{}{d}{B}\right)\sum _{l=1}^{\frac{n}{d}}\mu (l)\lfloor \frac{n}{dl}\rfloor \lfloor \frac{m}{dl}\rfloor$$
化到这步的时候我想的是
对前面整除分块，组合数列求和，内层套一个整除分块，内层套的时候用杜教筛 $\mu$ 的前缀和
算了一下复杂度，发现是 ${\sum }_{i=1}^{\sqrt{n}}\sqrt{\frac{n}{i}}=n^{3/4}$，加上杜教筛的 $O(n^{2/3})$，亲测大样例 $10s$
我们按照套路对上式继续化简，枚举 $dl$
$$\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor \sum _{l|T}\mu (l)\left(\genfrac{}{}{0ex}{}{T/l}{B}\right)$$
考虑直接筛后面一坨的前缀和就可以整除分块
令 $f(T)={\sum }_{l|T}\mu (l)\left(\genfrac{}{}{0ex}{}{T/l}{B}\right)$
$S(n)={\sum }_{i=1}^{n}f(n)$
考虑到有 $\mu$，卷上一个 $I$
那么就有
$$I\ast S(n)=\sum _{i=1}^n\sum _{l|i}\left(\genfrac{}{}{0ex}{}{i/l}{B}\right)[l=1]=\sum _{i=1}^n\left(\genfrac{}{}{0ex}{}{i}{B}\right)=\left(\genfrac{}{}{0ex}{}{n+1}{B+1}\right)$$
所以
$$I\ast S(n)=S(n)+\sum _{i=2}^{n}S(\lfloor \frac{n}{i}\rfloor )=\left(\genfrac{}{}{0ex}{}{n+1}{B+1}\right)$$
杜教筛即可，预处理往大的预处理，由于带 $log$，到 $3e6$ 就已经差不多了
复杂度是整除分块的 $O(\sqrt{n})$ 和杜教筛的 $O(n^{2/3})$
然后组合数还要用一下 $lucas$

#include<bits/stdc++.h>
#define cs const
#include<tr1/unordered_map>
using namespace std;
using namespace tr1;
typedef long long ll;
ll read(){
	ll cnt = 0, f = 1; char ch = 0;
	while(!isdigit(ch)){ ch = getchar(); if(ch == '-') f = -1; }
	while(isdigit(ch)) cnt = cnt*10 + (ch-'0'), ch = getchar();
	return cnt * f;
}
cs int Mod = 9990017;
ll n, m, B;
int fac[Mod + 5], ifac[Mod + 5];
ll add(ll a, ll b){ return (a + b) >= Mod ? (a + b) % Mod : a + b; }
ll mul(ll a, ll b){ if(a >= Mod) a %= Mod; if(b >= Mod) b %= Mod; ll r = a * b; if(r >= Mod) r %= Mod; return r; }
ll ksm(ll a, ll b){ ll ans = 1; for(;b;b>>=1,a=mul(a,a)) if(b&1) ans = mul(ans, a); return ans; }
ll C(ll n, ll m){ 
	if(n < 0 || m < 0 || n < m) return 0;
	return mul(fac[n], mul(ifac[m], ifac[n-m]));
}
ll Lucas(ll n, ll m){ 
	if(n < Mod && m < Mod) return C(n, m);
	return mul(Lucas(n / Mod, m / Mod), C(n % Mod, m % Mod));
}
cs int K = 5e6;
int prim[K + 5]; bool isp[K + 5]; 
int tot, mu[K + 5], f[K + 5], g[K + 5];
ll calc(ll x){ return Lucas(x + 1, B + 1); }
unordered_map<ll, int> MP;
ll F(ll x){
	if(x <= K) return f[x];
	if(MP.count(x)) return MP[x];
	ll ans = Lucas(x + 1, B + 1); 
	for(ll l = 2, r; l <= x; l = r + 1){
		ll v = x / l; r = x / v;
		ans = add(ans, Mod - mul(r-l+1, F(v)));
	} return MP[x] = ans;
}
void prework(){
	fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
	for(int i = 2; i < Mod; i++) fac[i] = mul(fac[i-1], i);
	ifac[Mod-1] = ksm(fac[Mod-1], Mod-2);
	for(int i = Mod-2; i >= 2; i--) ifac[i] = mul(ifac[i+1], i+1);
	mu[1] = 1;
	for(int i = 2; i <= K; i++){
		if(!isp[i]) prim[++tot] = i, mu[i] = -1;
		for(int j = 1; j <= tot; j++){
			if(i * prim[j] > K) break;
			isp[i * prim[j]] = 1;
			if(i % prim[j] == 0) break;
			mu[i * prim[j]] = mu[i] * mu[prim[j]];
		}
	}
	for(int i = 1; i <= K; i++) g[i] = C(i, B);
	for(int i = 1; i <= K; i++){
		for(int j = i; j <= K; j += i){
			f[j] = add(f[j], mul(mu[i], g[j / i]));
		}
	}
	for(int i = 1; i <= K; i++) f[i] = add(f[i-1], f[i]);
}
int main(){
	n = read(), m = read(), B = read();
	if(n > m) swap(n, m);
	prework();
	ll ans = 0, las = 0;
	for(ll l = 1, r; l <= n; l = r + 1){
		ll v1 = n / l, v2 = m / l; r = min(n / v1, m / v2);
		ll nx = F(r);
		ans = add(ans, mul(add(nx, Mod - las), mul(v1, v2)));
		las = nx;
	} cout << ans; return 0;
}