---

$\sum_{i = 1}^n\sum_{j = 1}^ni * j * gcd(i,j)$i=1∑n​j=1∑n​i∗j∗gcd(i,j)$=\sum_{d = 1}^{n}\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j = 1}^{\lfloor\frac{n}{d}\rfloor} i * j * d^3 * [gcd(i,j) = 1]$=d=1∑n​i=1∑⌊dn​⌋​j=1∑⌊dn​⌋​i∗j∗d3∗[gcd(i,j)=1]$=\sum_{d = 1}^{n}\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j = 1}^{\lfloor\frac{n}{d}\rfloor} i * j * d^3 * \sum_{p | gcd(i,j)}\mu(p)$=d=1∑n​i=1∑⌊dn​⌋​j=1∑⌊dn​⌋​i∗j∗d3∗p∣gcd(i,j)∑​μ(p)$=\sum_{d = 1}^{n}\sum_{p =1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\sum_{i = 1}^{\lfloor\frac{n}{d * p}\rfloor}\sum_{j = 1}^{\lfloor\frac{n}{d * p}\rfloor} i * j * d^3*p^2$=d=1∑n​p=1∑⌊dn​⌋​μ(p)i=1∑⌊d∗pn​⌋​j=1∑⌊d∗pn​⌋​i∗j∗d3∗p2$=\sum_{T = 1}^{n}\sum_{d | T}\mu(\frac{T}{d})*d^3*(\frac{T}{d})^2*f({\lfloor\frac{n}{T}\rfloor})^2,f(n) = \sum_{i = 1}^ni$=T=1∑n​d∣T∑​μ(dT​)∗d3∗(dT​)2∗f(⌊Tn​⌋)2,f(n)=i=1∑n​i$=\sum_{T = 1}^{n}T^2\sum_{d | T}\mu(\frac{T}{d})*d*f({\lfloor\frac{n}{T}\rfloor})^2$=T=1∑n​T2d∣T∑​μ(dT​)∗d∗f(⌊Tn​⌋)2$=\sum_{T = 1}^{n}T^2\phi(T)*f({\lfloor\frac{n}{T}\rfloor})^2$=T=1∑n​T2ϕ(T)∗f(⌊Tn​⌋)2
令$f(T) = T ^2 * \phi(T)$f(T)=T2∗ϕ(T)
$n &lt;= 10^7$n<=107的话可以用线性筛预处理$f$f函数前缀和，然后分块。
n 出到了 $10^{11}$1011可以用杜教筛筛出$f$f函数的前缀和，然后分块

杜教筛的大致过程：选一个 积性函数 $g$g，构造 $h = g * f$h=g∗f（$*$∗代表狄利克雷卷积），可知 h 也是一个积性函数。$h(x) = \sum_{d | x} g(d) * f(\frac{x}{d})$h(x)=d∣x∑​g(d)∗f(dx​)$\sum_{i = 1}^nh(i) = \sum_{i = 1}^n\sum_{d | i} g(d) * f(\frac{i}{d})$i=1∑n​h(i)=i=1∑n​d∣i∑​g(d)∗f(di​)$=\sum_{d = 1}^n\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor} g(d) * f(i)$=d=1∑n​i=1∑⌊dn​⌋​g(d)∗f(i)$令S(n) = \sum_{i = 1}^nf(i)：\sum_{i = 1}^nh(i) =\sum_{d = 1}^n\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor} g(d) * f(i) =\sum_{d = 1}^ng(d)*S({\lfloor\frac{n}{d}\rfloor}) = g(1) * S(n) + \sum_{d = 2}^ng(d)*S({\lfloor\frac{n}{d}\rfloor})$令S(n)=i=1∑n​f(i)：i=1∑n​h(i)=d=1∑n​i=1∑⌊dn​⌋​g(d)∗f(i)=d=1∑n​g(d)∗S(⌊dn​⌋)=g(1)∗S(n)+d=2∑n​g(d)∗S(⌊dn​⌋)$移项：g(1) * S(n) = \sum_{i = 1}^nh(i) - \sum_{d = 2}^ng(d)*S({\lfloor\frac{n}{d}\rfloor})$移项：g(1)∗S(n)=i=1∑n​h(i)−d=2∑n​g(d)∗S(⌊dn​⌋)若$h(i)$h(i)的前缀和可以在O（1）时间求得，则$S(n)$S(n)可以通过分块加记忆化搜索来求，因此要选一个积性函数$f$f使得 $h = g * f$h=g∗f且 $h$h 的前缀和容易求。
$观察g * f：\sum_{d | x}g(d) * f(\frac{x}{d}) = \sum_{d | x}g(d) * ({\frac{x}{d}})^2*\phi(\frac{x}{d})$观察g∗f：d∣x∑​g(d)∗f(dx​)=d∣x∑​g(d)∗(dx​)2∗ϕ(dx​)
当 $g(i)$g(i) 取 $i^2$i2时(方便去掉分母里的 $d^2$d2) $g * f = \sum_{d | x}d ^ 2 * \phi(\frac{x}{d}) =x ^ 3，前缀和公式为:\frac{x^2*(x + 1)^2}{4}$g∗f=d∣x∑​d2∗ϕ(dx​)=x3，前缀和公式为:4x2∗(x+1)2​
$\sum_{i = 1}^n i ^2 = \frac{i * (i + 1) * (2*i + 1)}{6}$i=1∑n​i2=6i∗(i+1)∗(2∗i+1)​
于是就可以杜教筛了，注意记忆化，先用线性筛预处理$\sqrt n$n​的前缀和，记忆化可以用map 也可以 unordered_map 或自己的hash，由于洛谷unordered_map编译不通过，这里还是用map

---

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 6e6;
map<ll,ll> mp;
bool ispri[maxn + 100];
int pri[maxn + 100],phi[maxn + 100],mu[maxn + 100];
ll sum[maxn + 100];
ll n,mod;
ll inv2,inv4,inv6;
ll fpow(ll a,ll b) {
	ll r = 1;
	while(b) {
		if(b & 1) r = r * a % mod;
		a = a * a % mod;
		b >>= 1;
	}
	return r;
}
void sieve(int n) {
	ispri[0] = ispri[1] = false;
	pri[0] = 0;mu[1] = 1;phi[1] = 1;
	for(int i = 2; i <= n; i++) {
		if(!ispri[i]) pri[++pri[0]] = i,mu[i] = -1,phi[i] = i - 1;
		for(int j = 1; j <= pri[0] && i * pri[j] <= n; j++) {
			ispri[i * pri[j]] = true;
			if(i % pri[j] == 0) {
				phi[i * pri[j]] = phi[i] * pri[j];
				break;
			}
			phi[i * pri[j]] = phi[i] * (pri[j] - 1);
			mu[i * pri[j]] = -mu[i];
		}
	}
	sum[0] = 0;
	for(int i = 1; i <= n; i++)
		sum[i] = (sum[i - 1] + phi[i] * (1ll * i * i % mod) % mod) % mod;
}
ll cal0(ll x) {
	x %= mod;
	return x * (x + 1) % mod * inv2 % mod;
}
ll cal1(ll x) {
	x %= mod;
	return x * x % mod * (x + 1) % mod * (x + 1) % mod * inv4 % mod;
}
ll cal2(ll x) {
	x %= mod;
	return x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;
}
ll getsum(ll p) {
	if(p <= maxn) return sum[p];
	if(mp[p]) return mp[p];
	ll l,r;
	ll res = cal1(p);
	for(l = 2; l <= p; l = r + 1) {
		r = (p / (p / l));
		ll t = getsum(p / l);
		res = (res + mod - t * ((cal2(r) - cal2(l - 1) + mod) % mod) % mod) % mod;
	}
	return mp[p] = (res + mod) % mod;
}
int main() {
	scanf("%lld%lld",&mod,&n);
	sieve(6000000);
	inv2 = fpow(2,mod - 2);
	inv4 = fpow(4,mod - 2);
	inv6 = fpow(6,mod - 2);
	ll l,r;
	ll ans = 0;
	for(l = 1; l <= n; l = r + 1) {
		r = (n / (n / l));
		ll t = (getsum(r) - getsum(l - 1) + mod) % mod;
		ll x = cal1(n / l);
		ans = (ans + x * t % mod) % mod;
	}
	printf("%lld\n",(ans + mod) % mod);
	return 0;
}