解法
题意就是要求：
\[ ( \sum_{i = 1}^n \sum_{j = 1}^n ij gcd(i, j) ) \mod p \]
言简意赅的题面非常舒适。出现了\(gcd\)可以考虑往反演方面思考。先进行例行的套路：枚举\(gcd\)。在这里为了书写方便，我就把取模省去。
\[ \sum_{d = 1}^n d \sum_{i = 1}^n \sum_{j = 1}^n ij [gcd(i, j) = d] \]
然后在后面消掉\(d\)这个玩意，同时约去两个\(d\)就会产生\(d^2\)，乘到外面即可：
\[ \sum_{d = 1}^n d^3 \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} ij [gcd(i, j) = 1] \]
发现后面的式子可以直接反演，考虑利用经典性质代换：
\[ \sum_{d = 1}^n d^3 \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} i \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} j \sum_{h = 1}^n \mu(h) [h|i][h|j] \]
调换式子位置：
\[ \sum_{d = 1}^n d^3 \sum_{h = 1}^n \mu(h) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} i[h|i] \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} j[h|j] \]
除掉两个\(h\)出来：
\[ \sum_{d = 1}^n d^3 \sum_{h = 1}^n \mu(h) h^2 \sum_{i = 1}^{\lfloor \frac{n}{hd} \rfloor} i \sum_{j = 1}^{\lfloor \frac{n}{hd} \rfloor} j \]
后面两个好像可以用高中数列求和用\(O(1)\)的时间求出，设\(calc(n) = \sum_{i = 1}^n i = \frac{n(n + 1)}{2}\)：
\[ \sum_{d = 1}^n d^3 \sum_{h = 1}^n \mu(h) h^2 calc(\lfloor \frac{n}{hd} \rfloor)^2 \]
设\(T = hd\)并化简：
\[ \sum_{d = 1}^n d^3 \sum_{T = 1}^n \mu(\lfloor \frac{T}{d} \rfloor) (\lfloor \frac{T}{d} \rfloor)^2 calc(\lfloor \frac{n}{T} \rfloor)^2 \\ \sum_{T = 1}^n T^2 calc(\lfloor \frac{n}{T} \rfloor)^2 \sum_{d|T} d \mu(\lfloor \frac{T}{d} \rfloor) \]
交换\(T\)和\(d\)位置时，要注意\(T = hd\)这个条件，故\(d|T\)需要被满足。我们发现后面这个式子：
\[ \sum_{d|T} d \mu(\lfloor \frac{T}{d} \rfloor) \]
这个式子是 Dirichlet 卷积的形式，简写为：
\[ id * \mu \]
对 Dirichlet 卷积熟悉的读者一眼可以看出来这两个函数的卷积为\(\varphi\)。这个结论一行就可以证完，如果不知道\( \epsilon = \mu * 1 \)怎么来的读者可以阅读这里，以下是证明：
\[ \varphi * 1 = id = id * \epsilon = id * \mu * 1 = id * \mu = \varphi \]
所以直接代换掉：
\[ \sum_{T = 1}^n calc(\lfloor \frac{n}{T} \rfloor)^2 \varphi(T) T^2 \]
已经很优秀了，用杜教筛搞下\(\varphi(T) T^2\)前缀和然后整除分块即可。
代码
// P3768.cpp
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll MAX_N = 5e6 + 2000;
ll n, mod, phi[MAX_N], prime[MAX_N], tot, inv2, inv6;
bool vis[MAX_N];
unordered_map<ll, ll> ump;
ll getPreQuad(ll num)
{
num %= mod;
return num * (num + 1) % mod * inv2 % mod;
}
ll getPreThird(ll num)
{
num %= mod;
return num * (num + 1) % mod * (num + num + 1) % mod * inv6 % mod;
}
void phi_sieve()
{
phi[1] = 1;
for (ll i = 2; i < MAX_N; i++)
{
if (!vis[i])
prime[++tot] = i, phi[i] = i - 1;
for (ll j = 1; j <= tot && 1ll * prime[j] * i < MAX_N; j++)
{
vis[i * prime[j]] = true;
if (i % prime[j] == 0)
{
phi[i * prime[j]] = phi[i] * prime[j] % mod;
break;
}
phi[i * prime[j]] = phi[i] * phi[prime[j]] % mod;
}
}
for (ll i = 2; i < MAX_N; i++)
phi[i] = (phi[i] * i % mod * i % mod + phi[i - 1]) % mod;
}
ll quick_pow(ll bas, ll tim)
{
ll ans = 1;
while (tim)
{
if (tim & 1)
ans = ans * bas % mod;
bas = bas * bas % mod;
tim >>= 1;
}
return ans;
}
ll sieve(ll num)
{
if (num < MAX_N)
return phi[num];
if (ump[num])
return ump[num];
ll ans = getPreQuad(num);
ans = ans * ans % mod;
for (ll i = 2, gx; i <= num; i = gx + 1)
{
gx = num / (num / i);
ll fact = (getPreThird(gx) - getPreThird(i - 1) + mod) % mod;
ans = (ans - fact * sieve(num / i) % mod + mod) % mod;
}
return ump[num] = (ans + mod) % mod;
}
signed main()
{
scanf("%lld%lld", &mod, &n);
inv2 = quick_pow(2, mod - 2), inv6 = quick_pow(6, mod - 2);
phi_sieve();
ll ans = 0;
for (ll i = 1, gx; i <= n; i = gx + 1)
{
gx = n / (n / i);
ll fact = getPreQuad(n / i);
fact = fact * fact % mod;
ll bfact = (sieve(gx) - sieve(i - 1) + mod) % mod;
ans = (ans + bfact * fact % mod) % mod;
}
printf("%lld", (ans + mod) % mod);
return 0;
}