【LG3768】简单的数学题

【LG3768】简单的数学题

题面

求

\[(\sum_{i=1}^n\sum_{j=1}^nij\text{gcd}(i,j))\text{mod}p \]

其中\(n\leq 10^{10},5\times 10^8\leq p \leq 1.1*10^9\)。

题解

推柿子：

\[\sum_{i=1}^n\sum_{j=1}^nij\text{gcd}(i,j)\\ =\sum_{d=1}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ijd^2[\gcd(i,j)==1]\\ =\sum_{d=1}d^3\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ij[\gcd(i,j)==1]\\ \]

令

\[f(d')=\sum_{i=1}^{n'}\sum_{j=1}^{n'}ij[\gcd(i,j)==d']\\ g(d')=\sum_{d'\mid x}f(x) \]

则

\[g(d')=\sum_{d'\mid x}\sum_{i=1}^{n'}\sum_{j=1}^{n'}ij[\gcd(i,j)==x]\\ =\sum_{i=1}^{n'}\sum_{j=1}^{n'}ij[d'|\gcd(i,j)]\\ =\sum_{i=1}^{\lfloor\frac {n'}{d'}\rfloor}\sum_{j=1}^{\lfloor\frac {n'}{d'}\rfloor}ijd'^2[1|\gcd(i,j)]\\ =S(\lfloor\frac {n'}{d'}\rfloor)^2d'^2 \]

其中\(S(x)=\sum_{i=1}^x i\)

那么

\[f(d')=\sum_{d'\mid x}S(\lfloor\frac {n'}x\rfloor)^2x^2\mu (x)\\ f(1)=\sum_{x=1}^{n'}S(\lfloor\frac {n'}x\rfloor)^2x^2\mu (x) \]

代回去

\[\sum_{d=1}^nd^3\sum_{x=1}^{\lfloor\frac nd\rfloor}S(\lfloor\frac {n}{xd}\rfloor)^2x^2\mu (x) \]

令\(Q=xd\)，

\[\sum_{d=1}^nd^3\sum_{x=1}^{\lfloor\frac nd\rfloor}S(\lfloor\frac {n}{Q}\rfloor)^2x^2\mu (x)\\ =\sum_{Q=1}^nS(\lfloor\frac {n}{Q}\rfloor)\sum_{d\mid Q}d^3(\frac Qd)^2\mu (\frac Qd)\\ =\sum_{Q=1}^nS(\lfloor\frac {n}{Q}\rfloor)^2\sum_{d\mid Q}dQ^2\mu (\frac Qd)\\ =\sum_{Q=1}^nS(\lfloor\frac {n}{Q}\rfloor)^2Q^2\varphi (Q) \]

现在我们如何求\(Q^2\varphi (Q)\)的前缀和呢？

令\(f=\text{id}^{2}\cdot \varphi\)，令\(g=\text{id}^2\cdot 1\)，

那么\(f*g=\text{id}^3,g=\text{id}^2\)。

就做完了。

#include <iostream> 
#include <cstdio> 
#include <cstdlib> 
#include <cstring> 
#include <cmath> 
#include <algorithm>
#include <map> 
using namespace std; 
const int MAX = 1e7, MAX_N = 1e7 + 5; 
bool nprime[MAX_N]; 
int prime[MAX_N], cnt, phi[MAX_N], f[MAX_N];
long long N; 
int Mod; 
int fpow(int x, int y) { 
	int res = 1;
	while (y) {
		if (y & 1) res = 1ll * res * x % Mod; 
		x = 1ll * x * x % Mod;
		y >>= 1; 
	} 
	return res; 
} 
void sieve() { 
	phi[1] = 1; 
	for (int i = 2; i <= MAX; i++) { 
		if (!nprime[i]) prime[++cnt] = i, phi[i] = i - 1; 
		for (int j = 1; j <= cnt && prime[j] * i <= MAX; j++) { 
			nprime[prime[j] * i] = 1; 
			if (i % prime[j] == 0) { phi[i * prime[j]] = 1ll * phi[i] * prime[j] % Mod; break; } 
			phi[i * prime[j]] = 1ll * phi[i] * phi[prime[j]] % Mod; 
		} 
	} 
	for (int i = 1; i <= MAX; i++) f[i] = (f[i - 1] + 1ll * phi[i] * i % Mod * i % Mod) % Mod; 
}
int inv2, inv6; 
int sqr(int n) { return 1ll * n * n % Mod; } 
int S(long long n) { return n %= Mod, n % Mod * (n + 1) % Mod * inv2 % Mod; } 
int g(long long n) { return n %= Mod, n % Mod * (n + 1) % Mod * (2 * n + 1) % Mod * inv6 % Mod; } 
map<long long, int> mp; 
int get_f(long long n) { 
	if (n <= MAX) return f[n]; 
	if (mp.find(n) != mp.end()) return mp[n]; 
	int res = sqr(S(n)); 
	for (long long l = 2, r; l <= n; l = r + 1) { 
		r = n / (n / l);
		res = (res - 1ll * (g(r) - g(l - 1) + Mod) % Mod * get_f(n / l) % Mod + Mod) % Mod; 
	} 
	return mp[n] = res; 
} 
int main () { 
#ifndef ONLINE_JUDGE 
    freopen("cpp.in", "r", stdin); 
#endif 
	cin >> Mod >> N; 
	sieve(); 
	inv2 = fpow(2, Mod - 2), inv6 = fpow(6, Mod - 2); 
	int ans = 0; 
	for (long long l = 1, r; l <= N; l = r + 1) { 
		r = N / (N / l); 
		ans = (ans + 1ll * sqr(S(N / l)) * ((get_f(r) - get_f(l - 1) + Mod) % Mod) % Mod) % Mod; 
	}
	printf("%d\n", ans); 
    return 0; 
}