【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;
}