题目描述
题面简单,我打下:
∑
i
=
1
n
∑
j
=
1
n
i
j
gcd
(
i
,
j
)
m
o
d
p
\sum_{i=1}^{n}\sum_{j=1}^{n}ij\gcd(i,j)\mod{p}
i=1∑nj=1∑nijgcd(i,j)modp
n最大1e10,p为质数
思路分析
拿到这种题,肯定是先推式子,暴力跑不过去的。
∑
i
=
1
n
∑
j
=
1
n
i
j
gcd
(
i
,
j
)
m
o
d
p
\sum_{i=1}^{n}\sum_{j=1}^{n}ij\gcd(i,j)\mod{p}
∑i=1n∑j=1nijgcd(i,j)modp
=
∑
d
=
1
n
d
3
∑
i
=
1
⌊
n
d
⌋
∑
j
=
1
⌊
n
d
⌋
[
gcd
(
i
,
j
)
=
1
]
=\sum_{d=1}^{n}{d^3}\sum_{i=1}^{⌊{n\over d}⌋}\sum_{j=1}^{⌊{n\over d}⌋}[\gcd(i,j)=1]
=∑d=1nd3∑i=1⌊dn⌋∑j=1⌊dn⌋[gcd(i,j)=1](改变枚举顺序)
=
∑
d
=
1
n
d
3
∑
i
=
1
⌊
n
d
⌋
i
∑
j
=
1
⌊
n
d
⌋
j
∑
x
∣
gcd
(
i
,
j
)
μ
(
x
)
=\sum_{d=1}^{n}{d^3}\sum_{i=1}^{⌊{n\over d}⌋}i\sum_{j=1}^{⌊{n\over d}⌋}j\sum_{x|\gcd(i,j)}\mu(x)
=∑d=1nd3∑i=1⌊dn⌋i∑j=1⌊dn⌋j∑x∣gcd(i,j)μ(x)(反演结论)
=
∑
d
=
1
n
d
3
∑
x
=
1
⌊
n
d
⌋
x
2
μ
(
x
)
∑
i
=
1
⌊
n
x
d
⌋
i
∑
i
=
1
⌊
n
x
d
⌋
j
=\sum_{d=1}^{n}{d^3}\sum_{x=1}^{⌊{n\over d}⌋}x^2\mu(x)\sum_{i=1}^{{⌊{n\over {xd}}⌋}}i\sum_{i=1}^{{⌊{n\over {xd}}⌋}}j
=∑d=1nd3∑x=1⌊dn⌋x2μ(x)∑i=1⌊xdn⌋i∑i=1⌊xdn⌋j(改变枚举顺序)
令sum(i)= ( 1 + i ) i 2 {(1+i)i\over2} 2(1+i)i,则,原式化为
=
∑
d
=
1
n
d
3
∑
x
=
1
⌊
n
d
⌋
μ
(
x
)
x
2
s
u
m
2
(
⌊
n
x
d
⌋
)
=\sum_{d=1}^{n}{d^3}\sum_{x=1}^{⌊{n\over d}⌋}\mu(x)x^2sum^2({⌊{n\over {xd}}⌋})
=∑d=1nd3∑x=1⌊dn⌋μ(x)x2sum2(⌊xdn⌋)
令T=xd
=
∑
T
=
1
n
T
2
s
u
m
2
(
⌊
n
T
⌋
)
∑
d
∣
T
μ
(
d
)
T
d
=\sum_{T=1}^{n}T^2sum^2(⌊{n\over {T}}⌋)\sum_{d|T}\mu(d){T\over d}
=∑T=1nT2sum2(⌊Tn⌋)∑d∣Tμ(d)dT
到此为止光从式子就推不下去了,但是,联想莫比乌斯函数的性质
(
μ
∗
i
d
)
=
ϕ
(\mu*id)=\phi
(μ∗id)=ϕ
也就是说,莫比乌斯函数和单位函数的迪利克雷卷积为欧拉函数。因此,原式可写为
∑
T
=
1
n
s
u
m
2
(
⌊
n
T
⌋
)
T
2
ϕ
(
T
)
\sum_{T=1}^{n}sum^2(⌊{n\over {T}}⌋)T^2\phi(T)
T=1∑nsum2(⌊Tn⌋)T2ϕ(T)
前面的sum函数部分可以用分块整除,那么,就需要考虑求出后面两个因式的部分和。
由于函数f(T)=T和欧拉函数都为积性函数,则,后面两个因式的乘积也为积性函数,再看看数据范围,故考虑杜教筛。
拿出杜教筛的套路式子:
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}^ng(d)S(⌊{n\over d}⌋)(答案式)
g(1)S(n)=i=1∑nh(i)−d=2∑ng(d)S(⌊dn⌋)(答案式)
其中
h
(
i
)
=
(
f
∗
g
)
(
i
)
h(i)=(f*g)(i)
h(i)=(f∗g)(i)
S
(
n
)
=
∑
i
=
1
n
f
(
i
)
S(n)=\sum_{i=1}^{n}f(i)
S(n)=∑i=1nf(i)
若能得出上式就能递归求解了。
重点在于g(x)的选择,我们看下我们要处理的式子:
f
(
x
)
=
x
2
ϕ
(
x
)
f(x)=x^2\phi(x)
f(x)=x2ϕ(x),带入上式得到:
g
(
1
)
S
(
n
)
=
∑
i
=
1
n
∑
d
∣
i
g
(
i
)
f
(
i
d
)
−
∑
d
=
2
n
g
(
d
)
S
(
⌊
n
d
⌋
)
(
1
)
g(1)S(n)=\sum_{i=1}^n\sum_{d|i}g(i)f({i\over d})-\sum_{d=2}^ng(d)S(⌊{n\over d}⌋)(1)
g(1)S(n)=i=1∑nd∣i∑g(i)f(di)−d=2∑ng(d)S(⌊dn⌋)(1)
单独看左边的一项:
∑
i
=
1
n
∑
d
∣
i
g
(
d
)
f
(
i
d
)
\sum_{i=1}^n\sum_{d|i}g(d)f({i\over d})
∑i=1n∑d∣ig(d)f(di)
=
∑
i
=
1
n
∑
d
∣
i
g
(
d
)
(
i
d
)
2
ϕ
(
i
d
)
(
2
)
=\sum_{i=1}^n\sum_{d|i}g(d)({i\over d})^2\phi({i\over d})(2)
=∑i=1n∑d∣ig(d)(di)2ϕ(di)(2)
联想欧拉函数性质:
∑
d
∣
n
ϕ
(
d
)
=
n
\sum_{d|n}\phi(d)=n
∑d∣nϕ(d)=n
因此我们考虑凑出上面的式子,也就是取
g
(
x
)
=
x
2
g(x)=x^2
g(x)=x2
我们回代:
对于(2)式:
原式=
∑
i
=
1
n
i
2
∑
d
∣
i
ϕ
(
d
)
=
∑
i
=
1
n
i
3
\sum_{i=1}^{n}i^2\sum_{d|i}\phi(d)=\sum_{i=1}^{n}i^3
∑i=1ni2∑d∣iϕ(d)=∑i=1ni3
带入(1)式:
S
(
n
)
=
∑
i
=
1
n
i
3
−
∑
d
=
2
n
g
(
d
)
S
(
⌊
n
d
⌋
)
(
g
(
1
)
=
1
)
S(n)=\sum_{i=1}^ni^3-\sum_{d=2}^ng(d)S(⌊{n\over d}⌋)(g(1)=1)
S(n)=i=1∑ni3−d=2∑ng(d)S(⌊dn⌋)(g(1)=1)
根据立方前缀和公式:
∑
i
=
1
n
i
3
=
(
1
+
2
+
…
…
+
n
)
2
=
s
u
m
2
(
n
)
\sum_{i=1}^ni^3=(1+2+……+n)^2=sum^2(n)
i=1∑ni3=(1+2+……+n)2=sum2(n)
就能杜教筛求解上面的式子了。
S
(
⌊
n
d
⌋
)
S(⌊{n\over d}⌋)
S(⌊dn⌋)整除分块+递归求解+线性筛预处理解决,g(d)的部分和,由于
g
(
x
)
=
x
2
g(x)=x^2
g(x)=x2,可以通过平方前缀和公式:
∑
i
=
1
n
i
2
=
n
(
n
+
1
)
(
2
n
+
1
)
6
\sum_{i=1}^ni^2={n(n+1)(2n+1)\over 6}
i=1∑ni2=6n(n+1)(2n+1)
转化为前缀和之差。
然后回代到上面的答案式就行了。
有点复杂的一题Orz。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
constexpr int maxn = 2e6 + 10;
template<typename T>
void read(T&x){
x=0;
int f=1;
char ch=getchar();
while(!isdigit(ch)){
if(ch=='-')f*=-1;
ch=getchar();
}
while(isdigit(ch)){
x=x*10+ch-'0';
ch=getchar();
}
x*=f;
}
template <typename T>
void write(T x)
{
if (x < 0)
putchar('-'), x = -x;
if (x > 9)
write(x / 10);
putchar(x % 10 + '0');
}
//=================================================
ll n, mod;
bitset<maxn> vis;
unordered_map<ll, int> dp;
int pri[maxn], tot;
int phi[maxn];
int pre_Sum[maxn];
int inv2, inv6;
inline int ksm(int a, int n)
{
int res = 1;
while (n)
{
if (n & 1)
res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
n >>= 1;
}
return res % mod;
}
inline void init()
{
inv2 = ksm(2, mod - 2);
inv6 = 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 = 0; j < tot && i * pri[j] < maxn; j++)
{
vis[i * pri[j]] = 1;
if (i % pri[j] == 0)
{
phi[i * pri[j]] = (1ll * phi[i]) * pri[j] % mod;
break;
}
phi[i * pri[j]] = (1ll * phi[i]) * (pri[j] - 1) % mod;
}
}
for (int i = 1; i < maxn; i++)
{
pre_Sum[i] = (((1ll * (i) * (i)) % mod) * phi[i]) % mod;
}
for (int i = 1; i < maxn; i++)
{
pre_Sum[i] = (pre_Sum[i] + pre_Sum[i - 1]) % mod;
}
}
inline int sum1(ll n)
{
n %= mod;
return (n * (n + 1) % mod * inv2 % mod) % mod;
}
inline int sum3(ll n)
{
n %= mod;
return (n * (n + 1) % mod * (2 * n + 1) % mod * inv6) % mod;
}
int L;
inline int sum2(ll n)
{
//if(L==1957302059)cerr<<"ok"<<endl;
//if(L==1957302059)cerr<<(n<maxn)<<endl,cerr<<n<<endl;
if (n < maxn)
return pre_Sum[n];
//if(L==1957302059)cerr<<"ok"<<endl;
if (dp.count(n))
return dp[n];
int tt = sum1(n) % mod;
int res = 1ll * tt * tt % mod;
for (ll l = 2, r; l <= n; l = r + 1)
{
//if(L==1957302059)cerr<<"**"<<l<<endl;
ll t = n / l;
r = n / t;
int p1 = ((sum3(r) - sum3(l - 1) % mod) + mod) % mod;
int p2 = sum2(t) % mod;
res = (1ll * (res - 1ll*p1 * p2 % mod) % mod);
}
return dp[n] = (res % mod + mod) % mod;
}
int solve()
{
int res = 0;
for (ll l = 1, r; l <= n; l = r + 1)
{
ll t = n / l;
r = n / t;
//cerr<<l<<" "<<r<<endl;
int p1 = sum1(t) % mod;
//if(l==1957302059) cerr<<l<<endl;
//cerr<<l<<" "<<r<<endl;
p1 = 1ll * p1 * p1 % mod;
int p2 = ((sum2(r) - sum2(l - 1)) % mod + mod) % mod;
//if(l==1957302059) cerr<<l<<endl;
//cerr<<l<<" "<<r<<endl;
res = (res % mod + 1ll * p1 * p2 % mod) % mod;
}
return (res % mod + mod) % mod;
}
signed main()
{
//freopen("in.txt", "r", stdin);
read(mod), read(n);
init();
write(solve());
return 0;
}