洛谷传送门
题目描述
由于出题人懒得写背景了,题目还是简单一点好。
输入一个整数 n n n和一个整数 p p p,你需要求出 ( ∑ i = 1 n ∑ j = 1 n i j g c d ( i , j ) ) m o d p (\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))~mod~p (∑i=1n∑j=1nijgcd(i,j)) mod p,其中 g c d ( a , b ) gcd(a,b) gcd(a,b)表示 a a a与 b b b的最大公约数。
刚才题面打错了,已修改
输入输出格式
输入格式:
一行两个整数 p p p、 n n n。
输出格式:
一行一个整数 ( ∑ i = 1 n ∑ j = 1 n i j g c d ( i , j ) ) m o d p (\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))~mod~p (∑i=1n∑j=1nijgcd(i,j)) mod p。
输入输出样例
输入样例#1:
998244353 2000
输出样例#1:
883968974
说明
对于20%的数据, n ≤ 1000 n \leq 1000 n≤1000。
对于30%的数据, n ≤ 5000 n \leq 5000 n≤5000。
对于60%的数据, n ≤ 1 0 6 n \leq 10^6 n≤106,时限1s。
对于另外20%的数据, n ≤ 1 0 9 n \leq 10^9 n≤109,时限3s。
对于最后20%的数据, n ≤ 1 0 10 n \leq 10^{10} n≤1010,时限6s。
对于100%的数据, 5 × 1 0 8 ≤ p ≤ 1.1 × 1 0 9 5 \times 10^8 \leq p \leq 1.1 \times 10^9 5×108≤p≤1.1×109且 p p p为质数。
解题分析
大力推式子:
∑
i
=
1
n
∑
j
=
1
n
i
j
g
c
d
(
i
,
j
)
=
∑
d
=
1
n
d
∑
i
=
1
n
∑
j
=
1
n
i
j
[
g
c
d
(
i
,
j
)
=
d
]
=
∑
d
=
1
n
d
∑
i
=
1
⌊
n
d
⌋
∑
j
=
1
⌊
n
d
⌋
i
j
d
2
[
g
c
d
(
i
,
j
)
=
1
]
=
∑
d
=
1
n
d
3
∑
i
=
1
⌊
n
d
⌋
∑
j
=
1
⌊
n
d
⌋
i
j
∑
p
∣
g
c
d
(
i
,
j
)
μ
(
p
)
=
∑
d
=
1
n
d
3
∑
p
=
1
⌊
n
d
⌋
μ
(
p
)
p
2
∑
i
=
1
⌊
n
d
p
⌋
∑
j
=
1
⌊
n
d
p
⌋
i
j
\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j) \\ =\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=d] \\ =\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ijd^2[gcd(i,j)=1] \\ =\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum_{p|gcd(i,j)}\mu(p) \\ =\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)p^2\sum_{i=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dp}\rfloor}ij
i=1∑nj=1∑nijgcd(i,j)=d=1∑ndi=1∑nj=1∑nij[gcd(i,j)=d]=d=1∑ndi=1∑⌊dn⌋j=1∑⌊dn⌋ijd2[gcd(i,j)=1]=d=1∑nd3i=1∑⌊dn⌋j=1∑⌊dn⌋ijp∣gcd(i,j)∑μ(p)=d=1∑nd3p=1∑⌊dn⌋μ(p)p2i=1∑⌊dpn⌋j=1∑⌊dpn⌋ij
设
s
(
n
)
=
∑
i
=
1
n
i
s(n)=\sum_{i=1}^ni
s(n)=∑i=1ni, 并设
T
=
d
p
T=dp
T=dp, 那么:
∑
d
=
1
n
d
3
∑
p
=
1
⌊
n
d
⌋
μ
(
p
)
p
2
∑
i
=
1
⌊
n
d
p
⌋
∑
j
=
1
⌊
n
d
p
⌋
i
j
=
∑
T
=
1
n
s
(
n
T
)
2
T
2
∑
d
∣
T
d
μ
(
T
d
)
\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)p^2\sum_{i=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dp}\rfloor}ij \\ =\sum_{T=1}^{n}s(\frac{n}{T})^2T^2\sum_{d|T}d\mu(\frac{T}{d})
d=1∑nd3p=1∑⌊dn⌋μ(p)p2i=1∑⌊dpn⌋j=1∑⌊dpn⌋ij=T=1∑ns(Tn)2T2d∣T∑dμ(dT)
而我们知道
ϕ
(
i
)
=
μ
(
i
)
∗
i
d
(
i
)
\phi(i)=\mu(i)*id(i)
ϕ(i)=μ(i)∗id(i), 那么实际上这玩意就是:
∑
T
=
1
n
s
(
n
T
)
2
T
2
ϕ
(
T
)
\sum_{T=1}^{n}s(\frac{n}{T})^2T^2\phi(T)
T=1∑ns(Tn)2T2ϕ(T)
前面可以下底分块, 关键在于后面的
f
(
x
)
=
x
2
ϕ
(
x
)
f(x)=x^2\phi(x)
f(x)=x2ϕ(x)上面, 这玩意该怎么科学筛出前缀和。
考虑到
ϕ
(
i
)
∗
I
=
i
d
(
i
)
\phi(i)*I=id(i)
ϕ(i)∗I=id(i), 我们可以设
g
(
x
)
=
x
2
g(x)=x^2
g(x)=x2, 那么
h
(
x
)
=
g
(
x
)
∗
f
(
x
)
=
∑
d
∣
x
d
2
∗
ϕ
(
x
d
)
(
x
d
)
2
=
∑
d
∣
x
x
2
ϕ
(
x
)
=
x
3
h(x)=g(x)*f(x) \\=\sum_{d|x}d^2*\phi(\frac{x}{d})(\frac{x}{d})^2 \\=\sum_{d|x}x^2\phi(x) \\ =x^3
h(x)=g(x)∗f(x)=d∣x∑d2∗ϕ(dx)(dx)2=d∣x∑x2ϕ(x)=x3
那么根据这个套路的式子:
g
(
1
)
S
(
n
)
=
∑
i
=
1
n
h
(
i
)
−
∑
d
=
2
n
g
(
d
)
S
(
⌊
n
d
⌋
)
g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)
g(1)S(n)=i=1∑nh(i)−d=2∑ng(d)S(⌊dn⌋)
和
∑
i
=
1
n
h
(
i
)
=
(
x
+
1
)
2
x
2
4
,
∑
i
=
1
n
g
(
i
)
=
n
∗
(
n
+
1
)
∗
(
2
n
+
1
)
6
\sum_{i=1}^{n}h(i)=\frac{(x+1)^2x^2}{4},\sum_{i=1}^{n}g(i)=\frac{n*(n+1)*(2n+1)}{6}
∑i=1nh(i)=4(x+1)2x2,∑i=1ng(i)=6n∗(n+1)∗(2n+1), 就可以科学筛出来了。
代码如下:
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <cstring>
#include <tr1/unordered_map>
#include <algorithm>
#define R register
#define IN inline
#define W while
#define gc getchar()
#define MX 5005000
#define ll long long
int pcnt;
std::tr1::unordered_map <ll, ll> mp;
int phi[MX], pri[MX];
ll sum[MX], val[MX], MOD, inv2, inv4, inv6;
bool npr[MX];
void get()
{
phi[1] = 1, val[1] = 1;
R int i, j, tar;
for (i = 2; i <= 5e6; ++i)
{
if (!npr[i]) pri[++pcnt] = i, phi[i] = i - 1;
for (j = 1; j <= pcnt; ++j)
{
tar = i * pri[j];
if (tar > 5e6) break;
npr[tar] = true;
if (!(i % pri[j])) {phi[tar] = phi[i] * pri[j]; break;}
phi[tar] = phi[i] * phi[pri[j]];
}
}
for (i = 1; i <= 5e6; ++i) sum[i] = (sum[i - 1] + 1ll * i * i % MOD * phi[i] % MOD) % MOD;
}
IN ll f1(ll val)
{return val % MOD * val % MOD * (val + 1) % MOD * (val + 1) % MOD * inv4 % MOD;}
IN ll f2(ll val)
{return val % MOD * (val + 1) % MOD * (2 * val % MOD + 1) % MOD * inv6 % MOD;}
IN ll fpow(ll base, ll tim)
{
ll ret = 1;
W (tim)
{
if (tim & 1)ret = ret * base % MOD;
base = base * base % MOD, tim >>= 1;
}
return ret;
}
ll solve(R ll bd)
{
if (bd <= 5e6) return sum[bd];
if (mp.count(bd)) return mp[bd];
ll ret = f1(bd);
for (R ll lef = 2, rig; lef <= bd; lef = rig + 1)
{
rig = bd / (bd / lef);
ret -= ((f2(rig) - f2(lef - 1) + MOD) % MOD) * solve(bd / lef) % MOD;
}
return mp[bd] = (ret + MOD) % MOD;
}
int main(void)
{
ll n, ans = 0, bd, mul;
scanf("%lld%lld", &MOD, &n);
get();
inv2 = fpow(2, MOD - 2);
inv4 = fpow(4, MOD - 2);
inv6 = fpow(6, MOD - 2);
for (R ll lef = 1, rig; lef <= n; lef = rig + 1)
{
rig = n / (n / lef); bd = n / lef;
mul = (1 + bd) % MOD * bd % MOD * inv2 % MOD;
mul = mul * mul % MOD;
(ans += mul * ((solve(rig) - solve(lef - 1) + MOD) % MOD) % MOD) %= MOD;
}
printf("%lld", ans);
}