\(\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);

}

巴特西

LG P3768 简单的数学题

\(\text{Problem}\)

\(\text{Solution}\)

\(\text{Code}\)

最新文章

热门文章