Description：

求： $$ \sum_{i=1}^n\sum_{j=1}^m\mu(i)\mu(j)\sigma_1(i,j) $$ $1\leqslant n,m\leqslant 500000$

Solution：

$$ \begin{align} &\sum_{i=1}^n\sum_{j=1}^m\mu(i)\mu(j)\sigma_1(i,j)\\ =&\sum_{i=1}^n\sum_{j=1}^m\mu(i)\mu(j)\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]x\frac jy\\ =&\sum_{i=1}^n\mu(i)\sum_{x|i}x\sum_{j=1}^m\mu(j)\sum_{y|j}[\gcd(x,y)=1]\frac jy\\ 设&f(x)=\sum_{j=1}^m\mu(j)\sum_{y|j}[\gcd(x,y)=1]\frac jy\\ 原式&=\sum_{i=1}^n\mu(i)\sum_{x|i}xf(x)可以O(n\ln n)求，所以我们的问题变成了求f\\ f(x)&=\sum_{j=1}^m\mu(j)\sum_{y|j}[\gcd(x,y)=1]\frac jy\\ &=\sum_{y=1}^m[\gcd(x,y)=1]\sum_{j=1}^m\mu(j)\frac jy\\ &=\sum_{y=1}^m[\gcd(x,y)=1]\sum_{j=1}^{\lfloor\frac m y\rfloor}\mu(yj)j\\ &然后根据\mu(i,j)=\mu(i)\times \mu(j)\times [\gcd(i,j)=1]\\ f(x)&=\sum_{y=1}^m[\gcd(x,y)=1]\mu(y)\sum_{j=1}^{\lfloor\frac my\rfloor}\mu(j)j[\gcd(j,y)=1]\\ \end{align} $$

设： $$ h(y)=\mu(y)\sum_{j=1}^{\lfloor\frac my\rfloor}\mu(j)j[\gcd(j,y)=1] $$ 那么： $$ \begin{align} f(x)&=\sum_{y=1}^m[\gcd(x,y)=1]h(y)\\ &=\sum_{y=1}^m\sum_{d|x.d|y}\mu(d)h(y)\\ &=\sum_{d|x}\mu(d)\sum_{y=1,d|y}^mh(y)\\ &=\sum_{d|x}\mu(d)\sum_{y=1}^{\lfloor\frac m d\rfloor}h(dy) \end{align} $$ 设： $$ l(x)=\sum_{y=1}^{\lfloor\frac mx\rfloor}h(xy) $$ 那么有： $$ f(x)=\sum_{d|x}\mu(d)l(d) $$ 那么如果可以求 $l$ ，就可以 $O(n\ln n)$ 算 $f$ ，于是我们的问题又变成了求 $l$ 。

如果我们知道了 $h$ ，那么 $l$ 可以暴力算，因为复杂度是调和级数。

于是我们的问题又变成了求 $h$ 。发现求 $h$ 实际上就是要求： $$ \begin{align} h(y)=&\sum_{j=1}^{\lfloor\frac m y\rfloor}\mu(j)j[\gcd(j,y)=1]\\ =&\sum_{j=1}^{\lfloor\frac m y\rfloor}\mu(j)j\sum_{d|j,d|y}\mu(d)\\ =&\sum_{d|y}\mu(d)\sum_{j=1,d|j}^{\lfloor\frac m y\rfloor}\mu(j)j\\ =&\sum_{d|y}\mu(d)\sum_{j=1}^{\lfloor\frac m {yd}\rfloor}\mu(dj)dj \end{align} $$ 再设： $$ p(y,d)=\sum_{j=1}^{\lfloor\frac m{yd}\rfloor}\mu(dj)dj $$ 这个我们可以预处理 $\mu(dj)dj$ 的前缀和求。

不过在不卡常的情况下 $h$ 其实可以 $O(n\log ^2n)$ 求，其中一个 $\log$ 用于求 $\gcd$ 常数极小。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,m;

#define MAXN 500010

#define N 500000

bool isprime[MAXN];

int prime[MAXN],tot = 0;

int mu[MAXN];

int gcd(int a,int b){return (b == 0 ? a : gcd(b,a % b));}

int res[MAXN];

int f[MAXN],l[MAXN],h[MAXN];

#define MOD 998244353

int calch(int y)

{

	int sum = 0;

	for(int j = 1;j <= m / y;++j)if(gcd(j,y) == 1)sum = (sum + 1ll * j * mu[j] % MOD * mu[y] % MOD) % MOD;

	return sum;

}

int main()

{

	for(int i = 2;i <= N;++i)isprime[i] = true;

	mu[1] = 1;

	for(int i = 2;i <= N;++i)

	{

		if(isprime[i])prime[++tot] = i,mu[i] = MOD - 1;

		for(int j = 1;j <= tot && i * prime[j] <= N;++j)

		{

			int k = i * prime[j];

			isprime[i * prime[j]] = false;

			if(i % prime[j] == 0){mu[k] = 0;break;}

			else mu[k] = (MOD - mu[i]) % MOD;

		}

	}

	scanf("%d%d",&n,&m);

	int ans = 0;

	for(int i = 1;i <= m;++i)h[i] = calch(i);

	for(int d = 1;d <= n;++d)for(int y = 1;y <= m / d;++y)l[d] = (l[d] + 1ll * h[y * d] * mu[d] % MOD) % MOD;

	for(int x = 1;x <= n;++x)for(int i = x;i <= n;i += x)f[i] = (f[i] + 1ll * l[x] * i % MOD) % MOD;

	for(int x = 1;x <= n;++x)for(int i = x;i <= n;i += x)res[i] = (res[i] + 1ll * f[x] * mu[i] % MOD) % MOD;

	for(int i = 1;i <= n;++i)ans = (ans + res[i]) % MOD;

	cout << ans << endl;

	return 0;

}