首先为了简洁 \(n = N,m = M\)。
然后 \(\mu(i)\mu(j)\sigma(ij)\) 这个式子，当 \(i\) 和 \(j\) 都没有平方质因子时， \[\begin{align*} & \mu(i)\mu(j)\sigma(ij) \\ = & \mu(i)\mu(j)\sigma(\dfrac i{\gcd(i,j)})\sigma(\dfrac j{\gcd(i,j)})\sigma(\text{gcd}^2(i,j)) \\ = & \mu(i)\mu(j)\sigma(i)\sigma(j)\dfrac{\sigma(\gcd^2(i,j))}{\sigma^2(\gcd(i,j))} \\ \end{align*}\]
当 \(i\) 或 \(j\) 有平方质因子时，显然 \(\mu(i)\mu(j) = 0\)，所以没有贡献。
于是 \[\begin{align*}
& \sum\limits_{i = 1}^n\sum\limits_{j = 1}^m\mu(i)\mu(j)\sigma(i)\sigma(j)\dfrac{\sigma(\gcd^2(i,j))}{\sigma^2(\gcd(i,j))} \\
= & \sum\limits_{d = 1}^{\min(n,m)}\dfrac{\sigma(\gcd^2(i,j))}{\sigma^2(\gcd(i,j))} \\
\cdot & \sum\limits_{id\le n}\sum\limits_{jd\le m}[\gcd(i,j)=1]\mu(id)\mu(jd)\sigma(id)\sigma(jd) \\
= & \sum\limits_{d = 1}^{\min(n,m)}\dfrac{\sigma(\gcd^2(i,j))}{\sigma^2(\gcd(i,j))} \\
\cdot & \sum\limits_{id\le n}\sum\limits_{jd\le m}\sum\limits_{k|\gcd(i,j)}\mu(k)\mu(id)\mu(jd)\sigma(id)\sigma(jd) \\
= & \sum\limits_{d = 1}^{\min(n,m)}\dfrac{\sigma(\gcd^2(i,j))}{\sigma^2(\gcd(i,j))} \\
\cdot & \sum\limits_{dk\le\min(n,m)}\mu(k)\sum\limits_{idk\le n}\sum\limits_{jdk\le m}\mu(idk)\mu(jdk)\sigma(idk)\sigma(jdk) \\
= & \sum\limits_{T = 1}^{\min(n,m)}f(n,T)f(m,T)g(T)
\end{align*}\]
其中 \(f(n,m) = \sum\limits_{im\le n}\mu(im)\sigma(im),g(n) = \sum\limits_{d|n}\dfrac{\sigma(d^2)}{\sigma^2(d)}\mu(\dfrac nd)\)。
求出所有的 \(f(n,T),f(m,T),g(T)\) 都是 \(O(n \log n)\) 的（认为 \(n,m\) 同阶），于是线性筛出 \(\mu(i),\sigma(i),\sigma(i^2)\) 这些东西就可以做了。
然鹅我处处求逆元常数过大（又没法递推）……于是要吸氧才能过……
代码： 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
using namespace std;
const int N = 5e5;
const long long mod = 998244353;
inline long long fpow(long long a,long long b)
{
long long ret = 1;
for(;b;b >>= 1)
(b & 1) && (ret = ret * a % mod),a = a * a % mod;
return ret;
}
int n,m;
int vis[N + 5],cnt,prime[N + 5];
long long sigma[N + 5],mn[N + 5],sqr[N + 5],sqrmn[N + 5],mu[N + 5];
long long f[2][N + 5],g[N + 5],t[N + 5];
long long ans;
int main()
{
mu[1] = sigma[1] = mn[1] = sqr[1] = sqrmn[1] = 1;
for(register int i = 2;i <= N;++i)
{
if(!vis[i])
mu[prime[++cnt] = i] = mod - 1,sigma[i] = 1 + i,mn[i] = 1 + i,sqr[i] = (1 + i + (long long)i * i % mod) % mod,sqrmn[i] = (1 + i + (long long)i * i % mod) % mod;
for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
{
sigma[i * prime[j]] = sigma[i] * fpow(mn[i],mod - 2) % mod * (mn[i] * prime[j] % mod + 1) % mod,mn[i * prime[j]] = (mn[i] * prime[j] % mod + 1) % mod;
sqr[i * prime[j]] = sqr[i] * fpow(sqrmn[i],mod - 2) % mod * ((sqrmn[i] * prime[j] % mod + 1) % mod * prime[j] % mod + 1) % mod,sqrmn[i * prime[j]] = ((sqrmn[i] * prime[j] % mod + 1) % mod * prime[j] % mod + 1) % mod;
break;
}
else
{
sigma[i * prime[j]] = sigma[i] * (prime[j] + 1) % mod,mn[i * prime[j]] = prime[j] + 1;
sqr[i * prime[j]] = sqr[i] * (((long long)prime[j] * prime[j] + prime[j] + 1) % mod) % mod,sqrmn[i * prime[j]] = ((long long)prime[j] * prime[j] % mod + prime[j] + 1) % mod;
mu[i * prime[j]] = mod - mu[i];
}
}
}
scanf("%d%d",&n,&m);
if(n > m)
swap(n,m);
for(register int i = 1;i <= n;++i)
for(register int j = i,k = 1;j <= n;j += i,++k)
g[j] = (g[j] + sqr[i] * fpow(sigma[i] * sigma[i] % mod,mod - 2) % mod * mu[k] % mod) % mod;
for(register int i = 1;i <= n;++i)
{
for(register int j = i;j <= n;j += i)
f[0][i] = (f[0][i] + mu[j] * sigma[j] % mod) % mod;
for(register int j = i;j <= m;j += i)
f[1][i] = (f[1][i] + mu[j] * sigma[j] % mod) % mod;
ans = (ans + f[0][i] * f[1][i] % mod * g[i] % mod) % mod;
}
printf("%lld\n",ans);
}