LG P3768 简单的数学题

\(\text{Problem}\)

求

\[\left(\sum_{i=1}^n \sum_{j=1}^n i j \gcd(i,j)\right) \bmod p \]

\(n \le 10^{10}，5 \times 10^8 \le p \le 1.1 \times 10^9\) 且 \(p \in \mathbb{P}\)

\(\text{Solution}\)

显然走欧拉反演

\[\begin{aligned} \sum_{i=1}^n \sum_{j=1}^n i j \gcd(i,j) &= \sum_{i=1}^n \sum_{j=1}^n i j \sum_{d|\gcd(i,j)} \varphi(d) \\ &= \sum_{d=1}^n \varphi(d) \sum_{d|i} i \sum_{d|j} j \\ &= \sum_{d=1}^n \varphi(d) \sum_{i=1}^{\lfloor \frac n d \rfloor} i \sum_{j=1}^{\lfloor \frac n d \rfloor} j \\ &= \sum_{d=1}^n d^2 \varphi(d) S^3 (\lfloor \frac n d \rfloor) \end{aligned} \]

\(\varphi\) 后面的部分可以用等差数列求和公式的平方得到（它恰恰连续自然数三次方和的求和公式）
这个式子显然可以数论分块（非常显然）
那么重点就是求 \(S(n)=\sum_{d=1}^n d^2 \varphi(d)\)
杜教筛即可

即考虑卷积 \(f * g\)，记 \(f(n)=n^2 \varphi(n)\)，令 \(g = ID^2\)

\[(f * g)(n) = \sum_{d|n} f(d) g(\frac{n}{d}) = \sum_{d|n} d^2 \varphi(d) \left(\frac{n}{d}\right)^2 = n^2 \sum_{d|n} \varphi(d) = n^3 \]

那么

\[\begin{aligned} g(1)S(n)=\sum_{i=1}^n (f*g)(n) - \sum_{i=2}^n g(i)S(\lfloor \frac n i \rfloor) \\ S(n) = \sum_{i=1}^n i^3 - \sum_{i=2}^n i^2 S(\lfloor \frac n i \rfloor) \end{aligned} \]

仍然可以数论分块，利用平方和与立方和公式快速计算

\(\text{Code}\)

#include<cstdio>
#include<tr1/unordered_map>
#define LL long long
#define maxn 5000000
#define N 5000005
using namespace std;

int vis[N], phi[N], prime[N], totp;
LL P, n, inv6, sf[N];
tr1::unordered_map<LL, LL> SF;

inline LL fpow(LL x, LL y)
{
	LL res = 1;
	for(; y; y >>= 1)
	{
		if (y & 1) res = res * x % P;
		x = x * x % P;
	}
	return res;
}
inline LL S2(LL n)
{
	n %= P;
	return n * (n + 1) % P * (n * 2 + 1) % P * inv6 % P;
}
inline LL S3(LL n)
{
	n %= P;
	return n * (n + 1) / 2 % P * (n * (n + 1) / 2 % P) % P;
}

inline void sieve()
{
	vis[1] = phi[1] = 1;
	for(register int i = 2; i <= maxn; i++)
	{
		if (!vis[i]) prime[++totp] = i, phi[i] = i - 1;
		for(register int j = 1; j <= totp && prime[j] * i <= maxn; j++)
		{
			vis[i * prime[j]] = 1;
			if (i % prime[j]) phi[i * prime[j]] = phi[i] * phi[prime[j]];
			else{phi[i * prime[j]] = phi[i] * prime[j]; break;}
		}
	}
	for(register int i = 1; i <= maxn; i++) sf[i] = (sf[i - 1] + (LL)i * i % P * phi[i] % P) % P;
}

LL SumF(LL n)
{
	if (n <= maxn) return sf[n];
	if (SF[n]) return SF[n];
	LL res = S3(n), r;
	for(register LL l = 2; l <= n; l = r + 1)
	{
		r = n / (n / l);
		res = (res - (S2(r) - S2(l - 1) + P) % P * SumF(n / l) % P + P) % P;
	}
	return SF[n] = res;
}

int main()
{
	scanf("%lld%lld", &P, &n);
	sieve();
	LL ans = 0, r; inv6 = fpow(6, P - 2);
	for(register LL l = 1; l <= n; l = r + 1)
	{
		r = n / (n / l);
		ans = (ans + S3(n / l) * (SumF(r) - SumF(l - 1) + P) % P) % P;
	}
	printf("%lld\n", ans);
}