[题解] 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)=nmprodlimits_{p mid nm} frac{p-1}{p}=nmprodlimits_{p mid n ext{或} p mid m} frac{p-1}{p})

右边也类似的，上面就是

(nmprodlimits_{pmid n} frac{p-1}{p}prodlimits_{p mid m} frac{p-1}{p} gcd(n,m))

下面

(gcd(n,m)prodlimits_{p mid gcd(n,m)} frac{p-1}{p}=gcd(n,m)prodlimits_{p mid n ext{且} pmid m} frac{p-1}{p})

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

---

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

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

套路的去枚举(T=kd)

[Ans=sumlimits_{T} left(sumlimits_{dmid T} frac{d}{varphi(d)}mu(frac{T}{d}) ight) left(sumlimits_{i=1}^{lfloorfrac{n}{T} floor}varphi(iT) ight)left(sumlimits_{j=1}^{lfloorfrac{m}{T} floor}varphi(jT) ight) ]

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

[Ans=sumlimits_{T} f(T)g(T,lfloor n/T floor)g(T,lfloor m/T floor) ]

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

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

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

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

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

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

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

(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 floor > 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
", solve(n,m));
  }
  return 0;
}