洛谷 3768

根据Crash的数字表格，很容易可以将式子化简为

[egin{aligned} Ans &= sumlimits_{i = 1}^n sumlimits_{j = 1} ij(i, j) \ &= sumlimits_{d = 1}^n d^3 sumlimits_{k = 1}^{leftlfloorfrac{n}{d} ight floor} mu(k) k^2 left( sumlimits_{i = 1}^{leftlfloorfrac{n}{kd} ight floor} i ight)^2 end{aligned} ]

感觉 (d, k) 放在一起式子无法继续化简，主要是有 (kd) 存在，故令 (T = kd) ，则有

[Ans = sumlimits_{T = 1}^n left( sumlimits_{i = 1}^{leftlfloorfrac{n}{T} ight floor} i ight)^2 T^2 sumlimits_{d | T} d mu(frac{T}{d}) ]

那么考虑整除分块，现在需要处理的是后半部分

通过观察（看题解）可以发现， (sumlimits_{d | T} d mu(frac{T}{d})) 可以看成 (mu * id) ，故可替换成 (phi(T)) ，则有

[Ans = sumlimits_{T = 1}^n left( sumlimits_{i = 1}^{leftlfloorfrac{n}{T} ight floor} i ight)^2 T^2 phi(T) ]

现在考虑将 (T^2 phi(T)) 部分用杜教筛解决

[h(n) = sumlimits_{d | n} d^2 phi(d) g(frac{n}{d}) ]

为了消除 (d^2) ，令 (g(frac{n}{d}) = n^2) ，则有

[h(n) = n^3 ]

故得

[S(n) = sumlimits_{i = 1}^n i^3 - sumlimits_{d = 2}^n d^2 S(leftlfloorfrac{n}{d} ight floor) ]

又（通过看题解）有一个知识点

[sumlimits_{i = 1}^n i^3 = left( sumlimits_{i = 1}^n i ight) ]

那么就可以直接杜教筛了

至于复杂度，将最外围的整除分块与杜教筛看为一体，故复杂度为 (O (n^{frac{2}{3}}))

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <tr1/unordered_map>

using namespace std;

typedef long long LL;

const int MAXN = 6e06 + 10;

LL MOD, N;

LL power (LL x, LL p) {
	LL cnt = 1;
	while (p) {
		if (p & 1)
			cnt = cnt * x % MOD;
		x = x * x % MOD;
		p >>= 1;
	}
	return cnt;
}
LL inv2, inv6;

int prime[MAXN];
int vis[MAXN]= {0};
int pcnt = 0;
LL phi[MAXN]= {0};
LL sumphi[MAXN]= {0};
int MAX = 6e06;
void linear_sieve () {
	phi[1] = 1;
	for (int i = 2; i <= MAX; i ++) {
		if (! vis[i]) {
			prime[++ pcnt] = i;
			phi[i] = i - 1;
		}
		for (int j = 1; j <= pcnt && i * prime[j] <= MAX; j ++) {
			vis[i * prime[j]] = 1;
			if (! (i % prime[j])) {
				phi[i * prime[j]] = phi[i] * 1ll * prime[j] % MOD;
				break;
			}
			phi[i * prime[j]] = phi[i] * 1ll * (prime[j] - 1) % MOD;
		}
	}
	for (int i = 1; i <= MAX; i ++)
		sumphi[i] = (sumphi[i - 1] + 1ll * i % MOD * 1ll * i % MOD * phi[i] % MOD) % MOD;
}

tr1::unordered_map<LL, LL> maphi;
inline LL sqr (LL x) {
	return x * x % MOD;
}
inline LL eqm (LL n) {
	return (n + 1) % MOD * (n % MOD) % MOD * (2 * n % MOD + 1) % MOD * inv6 % MOD;
}
inline LL oseqm (LL n) {
	return n % MOD * ((n + 1) % MOD) % MOD * inv2 % MOD;
}
LL phi_sieve (LL n) {
	if (n <= MAX)
		return sumphi[n];
	if (maphi[n])
		return maphi[n];
	LL total = sqr (oseqm (n));
	for (LL l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		total = (total - (eqm (r) - eqm (l - 1) + MOD) % MOD * phi_sieve (n / l) % MOD + MOD) % MOD;
	}
	return maphi[n] = total;
}

LL Solve () {
	LL ans = 0;
	for (LL l = 1, r; l <= N; l = r + 1) {
		r = N / (N / l);
		ans = (ans + sqr (oseqm (N / l)) * ((phi_sieve (r) - phi_sieve (l - 1) + MOD) % MOD) % MOD) % MOD;
	}
	return ans;
}

int main () {
	scanf ("%lld%lld", & MOD, & N);
	inv2 = power (2ll, MOD - 2), inv6 = power (6ll, MOD - 2);
	MAX = (int) min (1ll * MAX, N), linear_sieve ();
	LL ans = Solve ();
	cout << ans << endl;

	return 0;
}

/*
998244353 2000
*/

/*
1000000007 9786510294
*/