洛谷 P3768 ：简单的数学题（莫比乌斯反演 + 杜教筛）

---

$$\sum _{i=1}^n\sum _{j=1}^{n}i\ast j\ast gcd(i,j)$$$$=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }i\ast j\ast d^3\ast [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\ast j\ast d^3\ast \sum _{p|gcd(i,j)}\mu (p)$$$$=\sum _{d=1}^n\sum _{p=1}^{\lfloor \frac{n}{d}\rfloor }\mu (p)\sum _{i=1}^{\lfloor \frac{n}{d\ast p}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d\ast p}\rfloor }i\ast j\ast d^3\ast p^2$$$$=\sum _{T=1}^n\sum _{d|T}\mu (\frac{T}{d})\ast d^3\ast (\frac{T}{d}{)}^2\ast f(\lfloor \frac{n}{T}\rfloor {)}^2,f(n)=\sum _{i=1}^{n}i$$$$=\sum _{T=1}^n{T}^2\sum _{d|T}\mu (\frac{T}{d})\ast d\ast f(\lfloor \frac{n}{T}\rfloor {)}^2$$$$=\sum _{T=1}^n{T}^2\varphi (T)\ast f(\lfloor \frac{n}{T}\rfloor {)}^2$$
令$f(T)=T^2\ast \varphi (T)$
$n<=10^7$的话可以用线性筛预处理$f$函数前缀和，然后分块。
n 出到了 $10^{11}$可以用杜教筛筛出$f$函数的前缀和，然后分块

杜教筛的大致过程：选一个 积性函数 $g$，构造 $h=g\ast f$（$\ast$代表狄利克雷卷积），可知 h 也是一个积性函数。$$h(x)=\sum _{d|x}g(d)\ast f(\frac{x}{d})$$$$\sum _{i=1}^{n}h(i)=\sum _{i=1}^n\sum _{d|i}g(d)\ast f(\frac{i}{d})$$$$=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }g(d)\ast f(i)$$$$令S(n)=\sum _{i=1}^{n}f(i)：\sum _{i=1}^{n}h(i)=\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }g(d)\ast f(i)=\sum _{d=1}^{n}g(d)\ast S(\lfloor \frac{n}{d}\rfloor )=g(1)\ast S(n)+\sum _{d=2}^{n}g(d)\ast S(\lfloor \frac{n}{d}\rfloor )$$$$移项：g(1)\ast S(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)\ast S(\lfloor \frac{n}{d}\rfloor )$$若$h(i)$的前缀和可以在O（1）时间求得，则$S(n)$可以通过分块加记忆化搜索来求，因此要选一个积性函数$f$使得 $h=g\ast f$且 $h$ 的前缀和容易求。
$$观察g\ast f：\sum _{d|x}g(d)\ast f(\frac{x}{d})=\sum _{d|x}g(d)\ast (\frac{x}{d}{)}^2\ast \varphi (\frac{x}{d})$$
当 $g(i)$ 取 $i^2$时(方便去掉分母里的 $d^2$) $$g\ast f=\sum _{d|x}{d}^2\ast \varphi (\frac{x}{d})=x^3，前缀和公式为:\frac{x^2\ast (x+1{)}^2}{4}$$
$$\sum _{i=1}^n{i}^2=\frac{i\ast (i+1)\ast (2\ast i+1)}{6}$$
于是就可以杜教筛了，注意记忆化，先用线性筛预处理$\sqrt{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;
}