题意描述
求 (displaystyle(sum_{i=1}^{n}sum_{j=1}^{n} ij gcd(i,j) )mod p)
数据范围： (nleq 10^{10}, p 为质数)。
solution
推了好久才推出来的一道题。
先枚举一下 (gcd(i,j)) 可得：
(displaystylesum_{d=1}^{n} dsum_{i=1}^{n}sum_{j=1}^{n} ij[gcd(i,j) == d])
后面柿子可以提出来一个 (d) 变为：
(displaystyle sum_{d=1}^{n}d^3 sum_{i=1}^{{nover d}}sum_{j=1}^{nover d} ij[gcd(i,j) == 1])
莫比乌斯反演一下可以变为：
(displaystylesum_{ d=1}^{n} d^3 sum_{i=1}^{nover d}sum_{j=1}^{nover d}sum_{pmid gcd(i,j)} mu(p) ij)
(displaystylesum_{d=1}^{n} d^3sum_{p=1}^{nover d} mu(p)sum_{i=1}^{nover d}sum_{j=1}^{nover d} ij[pmid i,pmid j])
我们观察一下后面的 (displaystylesum_{i=1}^{nover d}sum_{j=1}^{nover d} ij[pmid i,pmid j]) 这个柿子，你会发现：
(i) 的取值可能为： ( p, 2p, 3p, 4p....{nover dp})
(j) 的取值可能为：( p,2p, 3p, 4p...{nover dp})
设 (S(n) = 1+2+3+4....+n), 那么 (displaystylesum_{i=1}^{nover d}sum_{j=1}^{nover d} ij[pmid i,pmid j] = p^2 S({nover dp})^2)。
带入原式可得：
(displaystylesum_{d=1}^{n} d^3 sum_{p=1}^{nover d} mu(p)p^2 S({nover dp})^2)
(dp) 化 (Q) 可得：
(displaystylesum_{Q=1}^{n} S({nover Q})^2 sum_{dmid Q} d^3 ({Qover d})^2 mu({Qover d}))
(= displaystyle sum_{Q=1}^{n} S({nover Q})^2 sum_{dmid Q} Q^2 d mu({Qover d}))
(= displaystylesum_{Q=1}^{n}S({nover Q})^2 Q^2 sum_{dmid Q} d mu({Qover d}))
然后你会发现 (displaystyle sum_{dmid Q} dmu({Qover d})) 他是个卷积的形式，由 (phi = mu * id) 可得： (displaystyle sum_{dmid Q} dmu({Qover d}) = phi(Q)), 代入原式可得：
(displaystylesum_{Q=1}^{n} S({nover Q}) ^2 Q^2 phi(Q))
前面那个可以拿整除分块来解决，设 (f(i) = i^2 phi(i)) ,现在关键是怎么求出来 (f(i)) 的前缀和。
推到这里自己就不会了，看了看题解总算懂了。
(n leq 10^{10}), 所以线性筛又萎了，在这个数据范围下，能做的好像只剩下杜教筛了。
套一下积性函数的柿子 ： (displaystyle h(n) = sum_{dmid n} f(d) g({nover d}))
因为 (f(i)) 里面有 (i^2) ，所以考虑构造 (g({nover i})) 把这一项消掉，那么就 设 (g(i) = i^2) 。
根据小学知识可得：(displaystyle sum_{i=1}^{n} g(i) = {n(n+1)(2n+1)over 2}), (g(i)) 的前缀和就可以 (O(1)) 的算出来了。
把 (g(i)) 代入原式可得：
(displaystyle h(n) = sum_{dmid n} d^2 phi(d) ({nover d})^2)
(displaystyle h(n) = n^2 sum_{dmid n} phi(d))
又因为 (displaystyle sum_{dmid n}phi(d) = n), 代入可得： (h(n) = n^3)
根据小学知识可得： $displaystylesum_{i=1}^n h(i) = ({n imes (n+1)over 2})^2 $, 这可也可以 (O(1)) 的来求。
然后使用杜教筛就可以在非线性的时间内把 (f(i)) 的前缀和算出来了。
最后的答案用数论分块计算一下即可，然后这道题就做完了QWQ。
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<map>
using namespace std;
#define int long long
using namespace std;
const int N = 1e7+10;
int p,n,tot,inv2,inv6,ans;
int prime[N],phi[N],f[N];
bool check[N];
map<int,int> vis;
inline int read()
{
int s = 0,w = 1; char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
return s * w;
}
int ksm(int a,int b)
{
int res = 1;
for(; b; b >>= 1)
{
if(b & 1) res = res * a % p;
a = a * a % p;
}
return res;
}
void YYCH()
{
phi[1] = 1;
for(int i = 2; i <= N-5; i++)
{
if(!check[i])
{
prime[++tot] = i;
phi[i] = i-1;
}
for(int j = 1; j <= tot && i * prime[j] <= N-5; j++)
{
check[i * prime[j]] = 1;
if(i % prime[j] == 0)
{
phi[i * prime[j]] = phi[i] * prime[j] % p;
break;
}
else phi[i * prime[j]] = phi[i] * phi[prime[j]] % p;
}
}
for(int i = 1; i <= N-5; i++) f[i] = i * i % p * phi[i] % p;
for(int i = 1; i <= N-5; i++) f[i] = (f[i] + f[i-1]) % p;
}
int sqr(int x)
{
x %= p;
return x * (x+1) % p * (2 * x % p + 1) % p * inv6 % p;
}
int Sum(int n)
{
if(n == 0) return 0;
n %= p;//这里 n*n 直接算的话会炸long long 所以要预先模一下
int tmp = n * (n + 1) % p * inv2 % p;
return tmp * tmp % p;
}
int get(int n)
{
if(n <= N-5) return f[n];
if(vis[n]) return vis[n];
int res = Sum(n);
for(int l = 2, r; l <= n; l = r+1)
{
r = min(n,n/(n/l));
res -= get(n/l) * ((sqr(r)-sqr(l-1) % p + p) % p) % p;
res = (res + p) % p;
}
vis[n] = (res % p + p) % p;
return vis[n];
}
signed main()
{
p = read(); n = read(); YYCH(); inv2 = ksm(2,p-2); inv6 = ksm(6,p-2);
for(int l = 1, r; l <= n; l = r+1)
{
r = min(n,n/(n/l));
ans = (ans + Sum(n/l) * ((get(r)-get(l-1) + p) % p) % p) % p;
}
printf("%lld
",ans);
return 0;
}