luogu P4240 毒瘤之神的考验

luogu P4240 毒瘤之神的考验

首先有一个重要的式子
$$?(ij)=\frac{?(i)?(j)\gcd (i,j)}{?(\gcd (i,j)}$$
这题就直接这么推
$=\sum _{d=1}\frac{d}{?(d)}\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}?(id)?(jd)[\gcd (i,j)=1]$
$=\sum _{d=1}\frac{d}{?(d)}\sum _{t=1}^{n/d}\mu (t)\sum _{i=1}^{n/dt}\sum _{j=1}^{m/dt}?(idt)?(jdt)$
设$T=dt$
$=\sum _{T=1}^n\sum _{d|T}\frac{\mu (d/t)d}{?(d)}\sum _{i=1}^{n/T}?(iT)\sum _{j=1}^{m/T}?(jT)$

把三部分分别设一下，可以得到
$=\sum _{T=1}^{n}f(T)g(T,n/T)g(T,m/T)$

$f$很容易预处理
考虑$g$,容易发现$g(T,n)=g(T,n?1)+?(Tn)$

然后再把整条式子设出来，

$h(a,b,n)=\sum _{T=1}^{n}f(T)g(T,a)g(T,b)$

考虑阈值分治，设一个阈值$B$
预处理$h(a,b,n),a\le B,b\le B$的答案

对于$n/T,m/T\ge B$的直接暴力跑
剩下的整除分块
时间复杂度
$O(n\ln n+n{B}^2+T(\sqrt{n}+\frac{n}{B})$
块大小我取的是$B=50$

code:

#include<bits/stdc++.h>
#define N 100005
#define B 50
#define mod 998244353
#define ll long long
using namespace std;
ll qpow(ll x, ll y) {
    ll ret = 1;for(; y; y >>= 1, x = x * x % mod) if(y & 1) ret = ret * x % mod;return ret;
}
int mu[N], phi[N], f[N], iphi[N];
vector<int> g[N], S[B + 5][B + 5]; 
int p[N], sz, vis[N];
void init(int n) {
    mu[1] = phi[1] = 1;for(int i = 2; i <= n; i ++) {
    if(!vis[i]) mu[i] = (mod - 1), p[++ sz] = i, phi[i] = i - 1;for(int j = 1; j <= sz && i * p[j] <= n; j ++) {
    vis[p[j] * i] = 1;if(i % p[j] == 0) {
    phi[i * p[j]] = 1ll * phi[i] * p[j] % mod;break;}phi[i * p[j]] = 1ll * phi[i] * phi[p[j]] % mod, mu[i * p[j]] = mod - mu[i];}}for(int i = 1; i <= n; i ++) iphi[i] = qpow(phi[i], mod - 2);for(int i = 1; i <= n; i ++)for(int j = 1; j <= n / i; j ++)f[i * j] = (f[i * j] + 1ll * i * mu[j] % mod * iphi[i] % mod) % mod;for(int T = 1; T <= n; T ++) {
    g[T].resize(n / T + 2);  g[T][0] = 0;for(int i = 1; i <= n / T; i ++) {
    g[T][i] = (g[T][i - 1] + phi[i * T]) % mod;}}for(int a = 1; a <= B; a ++)for(int b = 1; b <= B; b ++) {
     int lim = n / min(a, b);S[a][b].resize(lim + 2), S[a][b][0] = 0;for(int T = 1; T <= lim; T ++) {
    S[a][b][T] = (S[a][b][T - 1] + 1ll * f[T] * g[T][a] % mod * g[T][b] % mod) % mod;}}
}
ll calc(int n, int m) {
    ll ans = 0;for(int i = 1; i <= m / B; i ++) ans = (ans + 1ll * f[i] * g[i][n / i] % mod * g[i][m / i] % mod) % mod;for(int l = m / B + 1, r = 0; l <= n; l = r + 1) {
    r = min(n / (n / l), m / (m / l));ans = (ans + S[n / l][m / l][r] - S[n / l][m / l][l - 1] + mod) % mod;}return ans;
}
int n, m, t;
int main() {
    init(N - 5);scanf("%d", &t);while(t --) {
    scanf("%d%d", &n, &m); if(n > m) swap(n, m);printf("%lld\n", calc(n, m));}return 0;
}