( ext{Problem})
求
[left(sum_{i=1}^n sum_{j=1}^n i j gcd(i,j)
ight) mod p
]
(n le 10^{10}，5 imes 10^8 le p le 1.1 imes 10^9) 且 (p in mathbb{P})
( ext{Solution})
显然走欧拉反演
[egin{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
floor} i sum_{j=1}^{lfloor frac n d
floor} j \
&= sum_{d=1}^n d^2 varphi(d) S^3 (lfloor frac n d
floor)
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}
ight)^2 = n^2 sum_{d|n} varphi(d) = n^3
]
那么
[egin{aligned}
g(1)S(n)=sum_{i=1}^n (f*g)(n) - sum_{i=2}^n g(i)S(lfloor frac n i
floor) \
S(n) = sum_{i=1}^n i^3 - sum_{i=2}^n i^2 S(lfloor frac n i
floor)
end{aligned}
]
仍然可以数论分块，利用平方和与立方和公式快速计算
( ext{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
", ans);
}