Leaving aside this question, I found that after pushing more than a dozen times, I didn't even learn Dirichlet prefix and.

Dirichlet prefix and definition: \ (f*1=g \), we call \ (g \) the Dirichlet prefix and definition of \ (f \).

Given \ (f_1,f_2...f_n \), consider how to find the \ (g \) function in \ (O(n\log\log n) \).

Using the idea of high-dimensional prefix sum, we can regard each quality factor as one-dimensional, and then apply the solution of high-dimensional prefix sum.

For the binary prefix sum, we first accumulate one dimension, and then accumulate the other dimension. Figuratively, if we understand the binary prefix sum as a two-dimensional matrix, we first find the prefix sum for each row.

Then we can use the same idea to find the Dirichlet prefix and:

for (int i = 1; i <= cnt; ++ i)
  for (int j = 1; j * Prime[i] <= n; ++ j) f[j * Prime[i]] += f[j];

Dirichlet suffix and:

\(g_i = \ sum \ limits ^ {D \ Le n}{I \ vert d}f_d \), we call \ (g \) the Dirichlet suffix and of \ (f \).

Method: it is as like as two peas and de Lickley's prefix.

for (int i = 1; i <= cnt; ++ i)
  for (int j = n / Prime[i]; j; -- j) f[j] += f[j * Prime[i]];

Backward Dirichlet prefix and:

If we know \ (f*1=g \), how can we quickly find the \ (f \) function?

Solution: the difference array of prefix sum is the original array. Let's take two-dimensional prefix sum as an example. First, the rows are differentiated, and the resulting array is the prefix and sum of the corresponding columns. Just do the difference again.

But note that the difference should be done from large to small.

for (int i = 1; i <= cnt; ++ i)
  for (int j = n / Prime[i]; j; -- j) f[j * Prime[i]] -= f[j * Prime[i]];

Then we return to the question: please

\(1\le n,m\le 10^7\).

Change the original formula to

Then we found that this \ (\ sum \ limits {k = 1} ^ {ij} [ij\perp k] K \) is very difficult to move, but if \ (ij\perp k \), then \ (ij\perp ij-k \), so such \ (K \) always appear in pairs, so \ (\ sum \ limits {k = 1} ^ {ij} [ij\perp k] k = \ frac {\ varphi (ij) ij} {2} \), \ (ij=1 \) is an exception and requires special judgment.
How could I tell you I found it by typing my watch

So the original formula is (we ignore the special cases of \ (\ div 2 \) and \ (ij=1 \), because it doesn't matter)

It's still hard to move because \ (ij \) is big. At this time, only \ (\ varphi(ij)=\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))} \) can be used.

Then become

Make \ (f_1(i)=\sum\limits_{j=1}^{\lfloor\frac{n}{i}\rfloor}\varphi(ij)ij,f_2(i)=\sum\limits_{j=1}^{\lfloor\frac{m}{i}\rfloor}\varphi(ij)ij \), these two things are a Dirichlet suffix and can be preprocessed quickly.

Enumeration \ (dk \), converted to

Make \ (f (I) = \ frac {I ^ {n + 1} {\ varphi (I)}, G (I) = \ sum \ limits_ {D \ vert I} \ frac {d ^ {n + 1} {\ varphi (d)} \ mu (\ frac {I} {D}) \), obviously \ (g=f*\mu \).

Dirichlet convolution must be solved \ (O(n\log n) \), which cannot pass this problem. But \ (g=f*\mu \) means \ (f=g*1 \), which is to push back the Dirichlet prefix and. To sum up, the total time complexity \ (O(n\log\log n) \).

Then don't forget the last \ (+ 1 \) and then \ (\ div 2 \).

#include <cstdio>

const int mod = 998244353, N = 1e7 + 5;
int qpow(int a, int b) {
	int ret = 1ll;
	while (b) {
		if (b & 1) ret = 1ll * ret * a % mod;
		a = 1ll * a * a % mod, b >>= 1;
	}
	return ret;
}
int phi[N], inv[N], Prime[N / 10], cnt, f1[N], f2[N], g[N], n, m;
bool mark[10000005];
inline void add(int &x, const int y) {
	if ((x += y) >= mod) x -= mod;
}
inline void mns(int &x, const int y) {
	if ((x -= y) < 0) x += mod;
}

void init(int N) {
	g[1] = inv[1] = phi[1] = 1;
	for (int i = 2; i <= N; ++ i) {
		if (!mark[i]) Prime[++ cnt] = i, phi[i] = i - 1, g[i] = qpow(i, n + 1), inv[i] = qpow(i, mod - 2);
		for (int j = 1; i * Prime[j] <= N; ++ j) {
			mark[i * Prime[j]] = true, g[i * Prime[j]] = 1ll * g[i] * g[Prime[j]] % mod;
			inv[i * Prime[j]] = 1ll * inv[i] * inv[Prime[j]] % mod;
			if (i % Prime[j] == 0) {phi[i * Prime[j]] = phi[i] * Prime[j]; break;}
			phi[i * Prime[j]] = phi[i] * (Prime[j] - 1);
		}
	}
	for (int i = 1; i <= n; ++ i) g[i] = 1ll * g[i] * inv[phi[i]] % mod;
	for (int i = 1; i <= n; ++ i) f1[i] = 1ll * phi[i] * i % mod;
	for (int i = 1; i <= m; ++ i) f2[i] = 1ll * phi[i] * i % mod;
	for (register int i = cnt; i; -- i) 
		for (register int j = n / Prime[i]; j; -- j) mns(g[j * Prime[i]], g[j]);
	for (register int i = 1; i <= cnt && Prime[i] <= n; ++ i)
		for (register int j = n / Prime[i]; j; -- j) add(f1[j], f1[j * Prime[i]]);
	for (register int i = 1; i <= cnt && Prime[i] <= m; ++ i)
		for (register int j = m / Prime[i]; j; -- j) add(f2[j], f2[j * Prime[i]]);
}

int main() {
	scanf("%d%d", &n, &m);
	init(n > m ? n : m);
	int ans = 0;
	for (int i = 1; i <= n; ++ i) add(ans, 1ll * f1[i] * f2[i] % mod * g[i] % mod);
	printf("%d", ((int)((1ll * ans + 1) * qpow(2, mod - 2) % mod) + mod) % mod);
	return 0;
}