[题解] LuoguP4240 毒瘤之神的考验

https://www.luogu.com.cn/problem/P4240

---

毒瘤题......

中间那个\(\varphi\)就比较难搞......但注意到有这个柿子：

\[\varphi(nm)=\frac{\varphi(n)\varphi(m)\gcd(n,m)}{\varphi(\gcd(n,m))} \]

感性证一下，首先有

\(\varphi(nm)=nm\prod\limits_{p \mid nm} \frac{p-1}{p}=nm\prod\limits_{p \mid n \text{或} p \mid m} \frac{p-1}{p}\)

右边也类似的，上面就是

\(nm\prod\limits_{p\mid n} \frac{p-1}{p}\prod\limits_{p \mid m} \frac{p-1}{p} \gcd(n,m)\)

下面

\(\gcd(n,m)\prod\limits_{p \mid \gcd(n,m)} \frac{p-1}{p}=\gcd(n,m)\prod\limits_{p \mid n \text{且} p\mid m} \frac{p-1}{p}\)

最后一除\(\gcd(n,m)\)都没了，然后上下是个容斥的样子。

---

然后再来看题，就可以直接推柿子了。

\[\begin{aligned} Ans&=\sum\limits_{i=1}^n \sum\limits_{j=1}^m \frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))}\\&=\sum\limits_{d} \frac{d}{\varphi(d)} \sum\limits_{i=1}^{\lfloor n/d\rfloor} \sum\limits_{j=1}^{\lfloor m/d\rfloor}\varphi(id)\varphi(jd)[\gcd(i,j)=1]\\&=\sum\limits_{d} \frac{d}{\varphi(d)}\sum\limits_{i=1}^{\lfloor n/d\rfloor} \sum\limits_{j=1}^{\lfloor m/d\rfloor} \varphi(id)\varphi(jd) \sum\limits_{k \mid \gcd} \mu(k)\\&=\sum\limits_{d} \frac{d}{\varphi(d)}\sum\limits_{k} \mu(k) \sum\limits_{i=1}^{\lfloor\frac{n}{kd}\rfloor}\varphi(idk)\sum\limits_{j=1}^{\lfloor\frac{m}{kd}\rfloor}\varphi(jdk)\end{aligned} \]

套路的去枚举\(T=kd\)

\[Ans=\sum\limits_{T} \left(\sum\limits_{d\mid T} \frac{d}{\varphi(d)}\mu(\frac{T}{d})\right) \left(\sum\limits_{i=1}^{\lfloor\frac{n}{T}\rfloor}\varphi(iT)\right)\left(\sum\limits_{j=1}^{\lfloor\frac{m}{T}\rfloor}\varphi(jT)\right) \]

我们令\(f(n)=\sum\limits_{d \mid n} \frac{d}{\varphi(d)} \mu(n/d)\)，\(g(T,n)=\sum\limits_{i=1}^n \varphi(iT)\)，于是

\[Ans=\sum\limits_{T} f(T)g(T,\lfloor n/T\rfloor)g(T,\lfloor m/T\rfloor) \]

\(f\)很好办，筛出\(\varphi\)和\(\mu\)后调和级数复杂度预处理就好了。

\(g\)呢？里面有下取整，但还有一个参数\(T\)，咋办？

先不管别的，考虑把所有\(g\)都预处理出来，\(g(T,i)\)中有效的\(i\)只有\(\lfloor n/T\rfloor\)个

并且满足递推式\(g(T,i)=g(T,i-1)+\varphi(iT)\)。

现在我们求出了\(f\)和\(g\)，数论分块的希望还是很渺茫....

但我们还是可以往数论分块上想，考虑对\(\lfloor n/T\rfloor,\lfloor m/T\rfloor\)数论分块。

假设当前块为\([l,r]\)。我们考虑预处理一个东西\(h(a,b,n)=\sum\limits_{i=1}^n f(i)g(a,i)g(b,i)\)，那么我们就只要让答案加上\(h(r,\lfloor n/l\rfloor,\lfloor m/l\rfloor)-h(l-1,\lfloor n/l\rfloor,\lfloor m/l\rfloor)\)就可以了。

\(h\)可以非常显然的枚举\(a,b\)，然后有递推式\(h(a,b,i)=h(a,b,i-1)+f(i)g(i,a)g(i,b)\)。

预处理出所有\(h\)？复杂度原地爆炸...

同时注意到比较大的\(a,b\)的个数非常非常少，所以没有必要处理那么多，考虑预处理出\(h(1...B,1...B,n)\)。

那么对于大于\(B\)的\(a,b\)，暴力计算就好了。

至于\(B\)的取值，这里简单口胡一下复杂度，首先预处理是\(O(n \ln n + nB^2)\)的，然后单词询问，暴力算\(\lfloor n/i\rfloor > B\)的\(i\)对答案的贡献，这样的\(i\)差不多是\(n/B\)个，然后后面数论分块的复杂度就当\(O(\sqrt n)\)来算。

于是复杂度就是\(O(n \ln n + nB^2 + T(n/B + \sqrt{n}))\)，瞄一眼取\(B=n^{1/3}\)是会得到一个较为优秀的复杂度，大概是\(O(n \ln n + n^{5/3} + T(n^{2/3} + n^{1/2}))\)的。

#include <bits/stdc++.h>

using namespace std;

const int maxn = 1e5, N = maxn+10, mod = 998244353, maxb = 33;
typedef long long ll;

bool vis[N];
int ps[N], pn, mu[N], phi[N], inv[N];

int f[N];
vector<int> g[N], h[maxb + 10][maxb + 10];

void init() {
  int n = maxn; inv[1] = 1;
  for (int i = 2; i <= n; i++) inv[i] = 1ll * inv[mod%i] * (mod-mod/i) % mod;
  phi[1] = mu[1] = 1;
  for (int i = 2; i <= n; i++) {
    if (!vis[i]) {
      ps[pn++] = i;
      mu[i] = -1;
      phi[i] = i-1;
    }
    for (int j = 0; j < pn && i * ps[j] <= n; j++) {
      vis[i * ps[j]] = 1;
      if (i%ps[j] == 0) {
        mu[i*ps[j]] = 0;
        phi[i*ps[j]] = phi[i] * ps[j];
        break;
      }
      mu[i*ps[j]] = -mu[i];
      phi[i*ps[j]] = phi[i]*(ps[j]-1);
  }
  }
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j+=i) {
      f[j] += 1ll * mu[j/i] * i % mod * inv[phi[i]] % mod;
      if (f[j] >= mod) f[j] -= mod;
      if (f[j] < 0) f[j] += mod;
    }
  for (int i = 1; i <= n; i++) {
    int lim = n / i;
    g[i] = vector<int>(lim + 1);
    g[i][0] = 0;
    for (int j = 1; j <= lim; j++) {
      g[i][j] = g[i][j-1] + phi[i*j];
      if (g[i][j] >= mod) g[i][j] -= mod;
      if (g[i][j] < 0) g[i][j] += mod;
    }
  }
  
  for (int a = 1; a <= maxb; a++) {
    for (int b = 1; b <= maxb; b++) {
      int lim = min(n/a, n/b);
      auto &now = h[a][b];
      now = vector<int>(lim + 1);
      now[0] = 0;
      for (int i = 1; i <= lim; i++) {
        now[i] = now[i-1] + 1ll * g[i][a] * g[i][b] % mod * f[i] % mod;
        if (now[i] >= mod) now[i] -= mod;
        if (now[i] < 0) now[i] += mod;
      }
    }
  }
}

int solve(int n, int m) {
  int ans = 0, lst = 0; if (n > m) swap(n,m);
  for (int i = 1; ; i++) {
    if (n/i <= maxb && m/i <= maxb) { lst = i; break; }
    ans += 1ll * g[i][n/i] * g[i][m/i] % mod * f[i] % mod;
    if (ans >= mod) ans -= mod;
  }
  for (int l = lst, r = 0; l <= n; l = r+1) {
    r = min(n/(n/l), m/(m/l)); ans += (h[n/l][m/l][r] - h[n/l][m/l][l-1] + mod) % mod;
    if (ans >= mod) ans -= mod;
  }
  return ans;
}

int main() {
  init();
  int _, n, m; for (scanf("%d", &_); _; _--) {
    scanf("%d%d", &n, &m);
    printf("%d\n", solve(n,m));
  }
  return 0;
}