(color{#0066ff}{ 题目描述 })
由于出题人懒得写背景了，题目还是简单一点好。
输入一个整数n和一个整数p，你需要求出((sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p)，其中gcd(a,b)表示a与b的最大公约数。
(color{#0066ff}{输入格式})
一行两个整数p、n。
(color{#0066ff}{输出格式})
一行一个整数((sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p)。
(color{#0066ff}{输入样例})
998244353 2000
(color{#0066ff}{输出样例})
883968974
(color{#0066ff}{数据范围与提示})
对于20%的数据，(n leq 1000)。
对于30%的数据，(n leq 5000)。
对于60%的数据，(n leq 10^6)，时限1s。
对于另外20%的数据，(n leq 10^9)，时限3s。
对于最后20%的数据，(n leq 10^{10})，时限6s。
对于100%的数据，(5 imes 10^8 leq p leq 1.1 imes 10^9)且p为质数。
(color{#0066ff}{ 题解 })
这是一道神仙题qwq
开始推式子
题目要求
[sum_{i=1}^nsum_{j=1}^n i*j*gcd(i,j)
]
老规矩，枚举gcd
[sum_{d=1}^nsum_{i=1}^nsum_{j=1}^n [gcd(i,j)==d]i*j*d
]
把d提出来
[sum_{d=1}^n dsum_{i=1}^nsum_{j=1}^n [gcd(i,j)==d]i*j
]
后面把d弄出来
[sum_{d=1}^n d^3sum_{i=1}^{lfloor frac n d
floor}sum_{j=1}^{lfloor frac n d
floor} [gcd(i,j)==1]i*j
]
利用(mu * i = e)套进去
[sum_{d=1}^n d^3sum_{i=1}^{lfloor frac n d
floor}sum_{j=1}^{lfloor frac n d
floor} sum_{k|gcd(i,j)} mu(k)*i*j
]
即
[sum_{d=1}^n d^3sum_{k=1}^{lfloorfrac n d
floor}sum_{k|i} sum_{k|j} mu(k)*i*j
]
改为枚举倍数
[sum_{d=1}^n d^3sum_{k=1}^{lfloorfrac n d
floor} mu(k)*k^2sum _{i=1}^{lfloorfrac n {kd}
floor} sum _{j=1}^{lfloorfrac n {kd}
floor}i*j
]
看到这个式子，大为欣喜，于是
[令f(i)=mu(i) * i^2, g(i)=i^2
]
[h(n)=sum_{d|n} f(d)*g(frac n d)
]
[h(n)=n^2sum_{d|n} mu(d)=e(n)
]
于是以为这题切掉了？？？？？
我TM。。数据范围尼玛(10^9)
(how to O(sqrt n ^2))
显然要继续往下推
某巨佬有一手，pd换q的操作（原式为kd换q）
于是，先将原式变为
[sum_{d=1}^n sum_{k=1}^{lfloorfrac n d
floor} mu(k)*k^2 * d^3sum _{i=1}^{lfloorfrac n {kd}
floor} sum _{j=1}^{lfloorfrac n {kd}
floor}i*j
]
考虑枚举q，那么里面的kd就可以提出来了
[sum_{q=1}^n q^2sum_{k=1}^{lfloorfrac n d
floor} mu(k)* dsum _{i=1}^{lfloorfrac n {kd}
floor} sum _{j=1}^{lfloorfrac n {kd}
floor}i*j
]
这里面的d有点突兀qwq
[sum_{q=1}^n q^2sum_{k=1}^{lfloorfrac n d
floor} mu(k)* frac q k sum _{i=1}^{lfloorfrac n {q}
floor}i sum _{j=1}^{lfloorfrac n {q}
floor}j
]
上面的式子是有问题的，因为枚举的(q=kd)，那么k应该是q的因子，所以要改一下枚举，同时后面的可以写成平方形式
[sum_{q=1}^n q^2sum_{k|q} mu(k)* frac q k (sum _{i=1}^{lfloorfrac n {q}
floor}i)^2
]
然后可以把q提出来
[sum_{q=1}^n q^3sum_{k|q} frac {mu(k)} k (sum _{i=1}^{lfloorfrac n {q}
floor}i)^2
]
然后因为有
[sum_{k|n} frac {mu(k)} k =frac {varphi(n)} {n}
]
因此代回去
[sum_{q=1}^n q^3frac {varphi(q)} q (sum _{i=1}^{lfloorfrac n {q}
floor}i)^2
]
卧槽，继续
[sum_{q=1}^n varphi(q) * q^2 (sum _{i=1}^{lfloorfrac n {q}
floor}i)^2
]
我Fuck， (O(sqrt n)) 我就不信还过不了？？
[令f(i)=varphi(i) * i^2, g(i)=i^2
]
[h(n)=sum_{d|n} f(d)*g(frac n d)=sum_{d|n}varphi(d)*n^2=n^3 ，因为（varphi * i = id)
]
杜教筛+数列分块就行了
注意各种取模，少一个都WAqwq
#include<bits/stdc++.h>
#define int long long
#define LL long long
LL in() {
char ch; LL x = 0, f = 1;
while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
return x * f;
}
const int maxn = 4e6 + 10;
LL pri[maxn], phi[maxn], tot, mod;
LL six;
std::map<int, LL> mp;
bool vis[maxn];
LL ksm(LL x, LL y) {
LL re = 1LL;
while(y) {
if(y & 1) re = re * x % mod;
x = x * x % mod;
y >>= 1;
}
return re;
}
void predoit() {
mod = in();
six = ksm(6, mod - 2);
phi[1] = 1;
for(int i = 2; i < maxn; i++) {
if(!vis[i]) pri[++tot] = i, phi[i] = i - 1;
for(int j = 1; j <= tot && (LL)i * pri[j] < maxn; j++) {
vis[i * pri[j]] = true;
if(i % pri[j] == 0) {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
else phi[i * pri[j]] = phi[i] * (pri[j] - 1);
}
}
for(LL i = 2; i < maxn; i++) phi[i] = phi[i] * (i * i % mod) % mod;
for(LL i = 2; i < maxn; i++) (phi[i] += phi[i - 1]) %= mod;
}
LL getpfh(LL x) {
x %= mod;
LL ans = x;
ans = (ans * (x + 1)) % mod;
ans = (ans * ((x << 1LL) + 1)) % mod;
return ans * six % mod;
}
LL getlf(LL n) {
n %= mod;
n = (n * (n + 1) >> 1) % mod;
return n * n % mod;
}
LL getphi(LL n) {
if(n < maxn) return phi[n];
if(mp.count(n)) return mp[n];
LL ans = getlf(n);
for(LL l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
LL tot1 = ((getpfh(r) - getpfh(l - 1)) % mod + mod) % mod;
tot1 = tot1 * getphi(n / l) % mod;
ans = ((ans - tot1) % mod + mod) % mod;
}
return mp[n] = ans;
}
LL getsum(LL x) {
x %= mod;
LL ans = (x * (x + 1) >> 1LL) % mod;
return ans * ans % mod;
}
void work(LL n) {
LL ans = 0;
for(LL l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
LL tot1 = ((getphi(r) - getphi(l - 1)) % mod + mod) % mod;
tot1 = (tot1 * getsum(n / l)) % mod;
ans = (ans + tot1) % mod;
}
printf("%lld", (ans % mod + mod) % mod);
}
signed main() {
predoit();
work(in());
return 0;
}