文章目录
- 莫比乌斯函数
- 莫比乌斯反演
- 【CJOJ2512】gcd之和 -反演
- P2257 YY的GCD -反演
- 小清新数论 -杜教筛
- 求 ∑ i = 1 n ∑ j = 1 m l c m ( i , j ) 求\sum_{i=1}^n\sum_{j=1}^{m}lcm(i,j) 求∑i=1n∑j=1mlcm(i,j)
- 求 ∑ i = 1 n ∑ j = 1 m d ( i × j ) 求\sum_{i=1}^n\sum_{j=1}^{m}d(i\times j) 求∑i=1n∑j=1md(i×j)
- 小D的demo
- 2440: [中山市选2011]完全平方数
- HDU 6134
- BZOJ 3529 [Sdoi2014]数表
- 2019西安邀请赛B. Product
- 计蒜之道2019 1-D. 商汤AI园区的n个路口(困难)
- 插曲:二项式反演
- 插曲:集合反演
如有错误,忘大佬不吝赐教,及时指出。
积性函数好难啊,求大佬教我积性函数,或者推荐点博客。。
如果感觉有点卡, 请点这里。
莫比乌斯函数
参考:peng-ym
∑ d ∣ n μ ( d ) = [ n = 1 ] \sum_{d|n} \mu(d) = [n = 1] ∑d∣nμ(d)=[n=1] ( 拆 成 二 项 式 定 理 就 很 易 证 明 了 (拆成二项式定理就很易证明了 (拆成二项式定理就很易证明了
∑ d ∣ n ϕ ( d ) = n \sum_{d|n}\phi(d) = n ∑d∣nϕ(d)=n ( w i k i 上 的 证 明 挺 易 理 解 , 就 是 化 简 { 1 n , 2 n . . . n n } (wiki上的证明挺易理解,就是化简\{\frac 1n,\frac 2n...\frac nn\} (wiki上的证明挺易理解,就是化简{n1,n2...nn}
ϕ ( n ) = ∑ d ∣ n μ ( d ) n d \phi(n) = \sum_{d|n} \mu(d)\frac{n}{d} ϕ(n)=∑d∣nμ(d)dn ( 容 斥 过 程 ) (容斥过程) (容斥过程)
证明:
莫比乌斯反演
式 子 C : ∑ d ∣ g c d ( i , j ) μ ( d ) = [ g c d ( i , j ) = 1 ] 式子C: \sum_{d|gcd(i,j)} \mu(d) = [gcd(i,j) = 1] 式子C:∑d∣gcd(i,j)μ(d)=[gcd(i,j)=1]
式 子 B : F ( n ) = ∑ d ∣ n f ( d ) = > f ( n ) = ∑ d ∣ n u ( d ) F ( ⌊ n d ⌋ ) 式子B: F(n)=\sum_{d|n}f(d) => f(n)=\sum_{d|n}u(d)F(⌊\frac{n}{d}⌋) 式子B:F(n)=∑d∣nf(d)=>f(n)=∑d∣nu(d)F(⌊dn⌋)
式 子 A : F ( n ) = ∑ n ∣ d f ( d ) = > f ( n ) = ∑ n ∣ d u ( d n ) F ( d ) 式子A: F(n)=\sum_{n|d}f(d) => f(n)=\sum_{n|d}u(\frac{d}{n})F(d) 式子A:F(n)=∑n∣df(d)=>f(n)=∑n∣du(nd)F(d)
从莫比乌斯到欧拉:here
【CJOJ2512】gcd之和 -反演
求: a n s = ∑ i = 1 n ∑ j = 1 m g c d ( i , j ) n ≤ m ≤ 1 e 7 ans = \sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)\;\;\;n\leq m\leq1e7 ans=∑i=1n∑j=1mgcd(i,j)n≤m≤1e7
转换一下:
a
n
s
=
∑
d
=
1
n
d
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
d
]
=
∑
d
=
1
n
d
∑
i
=
1
n
d
∑
j
=
1
m
d
[
g
c
d
(
i
,
j
)
=
1
]
ans = \sum_{d=1}^{n}d \sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]=\sum_{d=1}^{n} d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[gcd(i,j)=1]
ans=d=1∑ndi=1∑nj=1∑m[gcd(i,j)=d]=d=1∑ndi=1∑dnj=1∑dm[gcd(i,j)=1]
法1:
f ( d ) = ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = d ] 并 且 F ( n ) = ∑ n ∣ d f ( d ) f(d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]并且F(n)=\sum_{n|d}f(d) f(d)=i=1∑nj=1∑m[gcd(i,j)=d]并且F(n)=n∣d∑f(d)
F ( n ) = ∑ i = 1 n ∑ j = 1 m [ d ∣ g c d ( i , j ) ] = ⌊ n d ⌋ × ⌊ m d ⌋ F(n)=\sum_{i=1}^{n}\sum_{j=1}^{m}[d|gcd(i,j)]=⌊\frac{n}{d}⌋\times ⌊ \frac{m}{d}⌋ F(n)=i=1∑nj=1∑m[d∣gcd(i,j)]=⌊dn⌋×⌊dm⌋
f ( x ) = ∑ x ∣ d m i n ( n , m ) u ( d x ) F ( d ) f(x)=\sum_{x|d}^{min(n,m)}u(\frac{d}{x})F(d) f(x)=x∣d∑min(n,m)u(xd)F(d)
f ( 1 ) = ∑ d = 1 m i n ( n , m ) u ( d ) F ( d ) = ∑ d = 1 m i n ( n , m ) u ( d ) ⌊ n d ⌋ × ⌊ m d ⌋ f(1)=\sum_{d=1}^{min(n,m)}u(d)F(d)=\sum_{d=1}^{min(n,m)}u(d)⌊\frac{n}{d}⌋\times ⌊ \frac{m}{d}⌋ f(1)=d=1∑min(n,m)u(d)F(d)=d=1∑min(n,m)u(d)⌊dn⌋×⌊dm⌋
a n s = ∑ d = 1 n d ∑ i = 1 n d u ( i ) n d × i m d × i ans = \sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}u(i)\frac{n}{d\times i}\frac{m}{d\times i} ans=d=1∑ndi=1∑dnu(i)d×ind×im
法2:(无敌易理解)
∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = 1 ] = ∑ i = 1 n ∑ j = 1 m ∑ d ∣ g c d ( i , j ) u ( d ) \sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=1]=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d|gcd(i,j)} u(d) i=1∑nj=1∑m[gcd(i,j)=1]=i=1∑nj=1∑md∣gcd(i,j)∑u(d)
∑ d = 1 n u ( d ) ∑ i = 1 n ∑ j = 1 m [ d ∣ g c d ( i , j ) ] = ∑ d = 1 n u ( d ) ⌊ n d ⌋ ⌊ m d ⌋ \sum_{d=1}^{n}u(d)\sum_{i=1}^{n}\sum_{j=1}^{m}[d|gcd(i,j)]=\sum_{d=1}^{n}u(d)⌊\frac{n}{d}⌋⌊\frac{m}{d}⌋ d=1∑nu(d)i=1∑nj=1∑m[d∣gcd(i,j)]=d=1∑nu(d)⌊dn⌋⌊dm⌋
a n s = ∑ d = 1 n d ∑ i = 1 n d u ( i ) n d × i m d × i ans = \sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}u(i)\frac{n}{d\times i}\frac{m}{d\times i} ans=d=1∑ndi=1∑dnu(i)d×ind×im
参考:这是一个神奇的Blog
AC_CODE: 和下面类似,就不写了。
P2257 YY的GCD -反演
题目:传送门
分析:
a
n
s
=
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
p
r
i
m
e
]
n
≤
m
≤
1
e
7
ans=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=prime]\;\;n\leq m\leq 1e7
ans=i=1∑nj=1∑m[gcd(i,j)=prime]n≤m≤1e7
把
p
r
i
m
e
prime
prime提到前面去:
a
n
s
=
∑
d
=
2
n
∑
i
=
1
n
d
∑
j
=
1
n
d
[
g
c
d
(
i
,
j
)
=
1
]
d
∈
p
r
i
m
e
ans=\sum_{d=2}^{n}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[gcd(i,j)=1]\;d\in prime
ans=d=2∑ni=1∑dnj=1∑dn[gcd(i,j)=1]d∈prime
莫比乌斯反演一下:
a
n
s
=
∑
d
=
2
n
∑
i
=
1
n
d
∑
j
=
1
n
d
∑
k
∣
g
c
d
μ
(
k
)
d
∈
p
r
i
m
e
ans=\sum_{d=2}^{n}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}\sum_{k|gcd}\mu(k)\;d\in prime
ans=d=2∑ni=1∑dnj=1∑dnk∣gcd∑μ(k)d∈prime
把
k
k
k提到前面来:
a
n
s
=
∑
d
=
2
n
∑
k
=
1
n
d
μ
(
k
)
n
k
d
m
k
d
d
∈
p
r
i
m
e
ans=\sum_{d=2}^{n}\sum_{k=1}^{\frac{n}{d}}\mu(k)\frac{n}{kd}\frac{m}{kd}\;d\in prime
ans=d=2∑nk=1∑dnμ(k)kdnkdmd∈prime
这样写会超时,因为多组数据,还得优化,就是交换求和顺序,令 T = k d : T=kd: T=kd:
a n s = ∑ T = 2 n n T m T ∑ d ∣ T μ ( T d ) d ∈ p r i m e ans=\sum_{T=2}^{n}\frac{n}{T}\frac{m}{T}\sum_{d|T}\mu(\frac{T}{d})\;d\in prime ans=T=2∑nTnTmd∣T∑μ(dT)d∈prime
O ( n l o g ( n ) ) O(nlog(n)) O(nlog(n))预处理这个式子 ∑ d ∣ T μ ( T d ) \sum_{d|T}\mu(\frac{T}{d}) ∑d∣Tμ(dT):枚举每个素数,给它的 x x x倍数加上 μ ( x ) \mu(x) μ(x),然后前缀和。
这个式子也可以
O
(
n
)
O(n)
O(n)预处理:
AC_CODE:
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MXN = 1e7+ 7;
const int mod = 998244353;
LL n, m;
int noprime[MXN], pp[MXN], pcnt;
int phi[MXN], mu[MXN];
LL pre_mu[MXN];
void init_prime() {
noprime[0] = noprime[1] = 1;
mu[1] = 1; phi[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, phi[i] = i-1, mu[i] = -1;
for(int j = 0; j < pcnt && pp[j] * i < MXN; ++j) {
noprime[pp[j]*i] = 1;
//phi[pp[j]*i] = (pp[j]-1)*phi[i];
mu[pp[j]*i] = -mu[i];
if(i % pp[j] == 0) {
//phi[pp[j]*i] = pp[j]*phi[i];
mu[pp[j]*i] = 0;
break;
}
}
}
}
void init_ans() {
for(int i = 0; i < pcnt; ++i) {
for(int j = 1; j*pp[i] < MXN; ++j) {
pre_mu[j*pp[i]] += mu[j];
}
}
for(int i = 2; i < MXN; ++i) pre_mu[i] += pre_mu[i-1];
}
int main() {
init_prime();
init_ans();
int tim; scanf("%d", &tim);
while(tim --) {
scanf("%lld%lld", &n, &m);
if(n > m) swap(n, m);
LL ans = 0, tmp;
for(LL L = 2, R; L <= n; L = R + 1) {
R = min(n/(n/L),m/(m/L));
tmp = pre_mu[R]-pre_mu[L-1];
ans = (ans + tmp*(n/L)*(m/L));
}
printf("%lld\n", ans);
}
return 0;
}
小清新数论 -杜教筛
题目:wannafly winter camp day3 F
a
n
s
=
∑
i
=
1
n
∑
j
=
1
n
μ
(
g
c
d
(
i
,
j
)
)
n
≤
1
e
11
ans=\sum_{i=1}^{n}\sum_{j=1}^{n}\mu(gcd(i,j))\;\;n\leq1e11
ans=i=1∑nj=1∑nμ(gcd(i,j))n≤1e11
这个式子与
g
c
d
gcd
gcd之和那题非常相似,如果数据范围是
1
e
7
1e7
1e7的话,这两题代码就基本一样了。可这题数据范围
1
e
11
1e11
1e11,就用杜教筛咯。
杜教筛学习博客:here
方法1:利用欧拉函数前缀和
a
n
s
=
∑
i
=
1
n
∑
j
=
1
n
μ
(
g
c
d
(
i
,
j
)
)
n
≤
1
e
11
ans=\sum_{i=1}^{n}\sum_{j=1}^{n}\mu(gcd(i,j))\;\;n\leq1e11
ans=i=1∑nj=1∑nμ(gcd(i,j))n≤1e11
a n s = ∑ d = 1 n μ ( d ) ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = d ] = ∑ d = 1 n μ ( d ) ∑ i = 1 n d ∑ j = 1 n d [ g c d ( i , j ) = 1 ] ans = \sum_{d=1}^{n}\mu(d) \sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]=\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[gcd(i,j)=1] ans=d=1∑nμ(d)i=1∑nj=1∑m[gcd(i,j)=d]=d=1∑nμ(d)i=1∑dnj=1∑dn[gcd(i,j)=1]
F ( n ) = ∑ i = 1 n ∑ j = 1 n [ g c d ( i , j ) = 1 ] , P ( n ) = ∑ i = 1 n ϕ ( i ) = ∑ i = 1 n ∑ j = 1 i [ g c d ( i , j ) = 1 ] F(n)=\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=1],\;P(n)=\sum_{i=1}^{n}\phi(i)=\sum_{i=1}^{n}\sum_{j=1}^{i}[gcd(i,j)=1] F(n)=i=1∑nj=1∑n[gcd(i,j)=1],P(n)=i=1∑nϕ(i)=i=1∑nj=1∑i[gcd(i,j)=1]
F ( n ) = P ( n ) × 2 − 1 F(n)=P(n)\times2 -1 F(n)=P(n)×2−1
a
n
s
=
∑
d
=
1
n
μ
(
d
)
×
(
P
(
n
d
)
×
2
−
1
)
ans = \sum_{d=1}^{n}\mu(d)\times(P(\frac{n}{d})\times 2-1)
ans=d=1∑nμ(d)×(P(dn)×2−1)、
方法2:利用莫比乌斯反演和迪利克雷卷积
a
n
s
=
∑
i
=
1
n
∑
j
=
1
n
μ
(
g
c
d
(
i
,
j
)
)
n
≤
1
e
11
ans=\sum_{i=1}^{n}\sum_{j=1}^{n}\mu(gcd(i,j))\;\;n\leq1e11
ans=i=1∑nj=1∑nμ(gcd(i,j))n≤1e11
化简一下:
a
n
s
=
∑
d
=
1
n
μ
(
d
)
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
d
]
=
∑
d
=
1
n
μ
(
d
)
∑
i
=
1
n
d
∑
j
=
1
n
d
[
g
c
d
(
i
,
j
)
=
1
]
ans = \sum_{d=1}^{n}\mu(d) \sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]=\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[gcd(i,j)=1]
ans=d=1∑nμ(d)i=1∑nj=1∑m[gcd(i,j)=d]=d=1∑nμ(d)i=1∑dnj=1∑dn[gcd(i,j)=1]
反演一下:
a
n
s
=
∑
d
=
1
n
μ
(
d
)
∑
i
=
1
n
d
∑
j
=
1
n
d
∑
k
∣
g
c
d
(
i
,
j
)
μ
(
k
)
ans=\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}\sum_{k|gcd(i,j)}\mu(k)
ans=d=1∑nμ(d)i=1∑dnj=1∑dnk∣gcd(i,j)∑μ(k)
把
k
k
k提出来:
a
n
s
=
∑
d
=
1
n
μ
(
d
)
∑
k
=
1
n
d
μ
(
k
)
(
n
k
d
)
2
ans=\sum_{d=1}^{n}\mu(d)\sum_{k=1}^{\frac{n}{d}}\mu(k)(\frac n{kd})^2
ans=d=1∑nμ(d)k=1∑dnμ(k)(kdn)2
又到了熟悉的步骤:令
T
=
k
d
T=kd
T=kd:
a
n
s
=
∑
T
=
1
n
(
n
T
)
2
∑
d
∣
T
μ
(
d
)
×
μ
(
T
d
)
=
∑
T
=
1
n
(
n
T
)
2
(
μ
∗
μ
)
(
T
)
ans=\sum_{T=1}^n(\frac nT)^2\sum_{d|T}\mu(d)\times \mu(\frac Td)=\sum_{T=1}^n(\frac nT)^2(\mu*\mu)(T)
ans=T=1∑n(Tn)2d∣T∑μ(d)×μ(dT)=T=1∑n(Tn)2(μ∗μ)(T)
接着问题变成求下面这个式子的前缀和:
∑
i
=
1
n
(
μ
∗
μ
)
(
i
)
\sum_{i=1}^n (\mu*\mu)(i)
i=1∑n(μ∗μ)(i)
AC_CODE1:
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MXN = 1e7+ 7;
const LL mod = 998244353;
const LL inf = 1e18;
LL n, m, inv2;
int noprime[MXN], pp[MXN], pcnt;
int phi[MXN], mu[MXN];
LL pre_mu[MXN], pre_phi[MXN];
//unordered_map<LL, LL> mp1, mp2;
void init_prime() {
noprime[0] = noprime[1] = 1;
mu[1] = 1; phi[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, phi[i] = i-1, mu[i] = -1;
for(int j = 0; j < pcnt && pp[j] * i < MXN; ++j) {
noprime[pp[j]*i] = 1;
phi[pp[j]*i] = (pp[j]-1)*phi[i];
mu[pp[j]*i] = -mu[i];
if(i % pp[j] == 0) {
phi[pp[j]*i] = pp[j]*phi[i];
mu[pp[j]*i] = 0;
break;
}
}
}
for(int i = 1; i < MXN; ++i) {
pre_mu[i] = pre_mu[i-1] + mu[i];
pre_phi[i] = (pre_phi[i-1] + phi[i])%mod;
}
}
struct Hash_map{
static const int mask=0x7fffff;
LL p[mask+1],q[mask+1];
void clear(){
memset(q,0,sizeof(q));
}
LL& operator [](LL k){
LL i;
for(i=k&mask;q[i]&&p[i]!=k;i=(i+1)&mask);
p[i]=k;
return q[i];
}
}mp1, mp2;
LL solve_mu(LL n) {
if(n < MXN) return pre_mu[n];
if(mp1[n] == inf) return 0;
if(mp1[n]) return mp1[n];
LL ans = 1;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans -= (R-L+1LL)%mod*solve_mu(n/L);
}
mp1[n] = ans;
if(mp1[n] == 0) mp1[n] = inf;
return ans;
}
LL solve_phi(LL n) {
if(n < MXN) return pre_phi[n];
if(mp2[n]) return mp2[n];
LL ans = n%mod*(n+1)%mod*inv2%mod;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans - (R-L+1LL)%mod*solve_phi(n/L)%mod)%mod;
}
ans = (ans + mod) % mod;
mp2[n] = ans;
return ans;
}
int main() {
inv2 = 499122177;
init_prime();
scanf("%lld", &n);
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (solve_mu(R)-solve_mu(L-1))%mod*(solve_phi(n/L)*2LL%mod-1LL)%mod+mod)%mod;
}
printf("%lld\n", (ans+mod)%mod);
return 0;
}
AC_CODE2:
#include<bits/stdc++.h>
#define fi first
#define se second
#define iis std::ios::sync_with_stdio(false);cin.tie(0)
#define pb push_back
#define o2(x) (x)*(x)
using namespace std;
typedef long long LL;
typedef pair<int, int> pii;
typedef pair<LL, LL> pll;
const int INF = 0x3f3f3f3f;
const LL MOD = 998244353;
const int MXN = 5e6 + 6;
LL n, q, m;
int noprime[MXN], pp[MXN/2], pcnt;
int mu[MXN], phi[MXN];
int pre_mu[MXN];//莫比乌斯函数的前缀和
LL mumu[MXN];//莫比乌斯函数卷积莫比乌斯函数的前缀和
unordered_map<LL,LL>mp1,mp2;
void init_rime() {
noprime[0] = noprime[1] = 1;
mu[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i]=-1;
for(int j = 0; j < pcnt && i*pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
mu[i*pp[j]] = -mu[i];
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
break;
}
}
}
for(int i = 1; i < MXN; ++i) pre_mu[i] = pre_mu[i-1] + mu[i], mumu[i] = mu[i];
for(int i = 2; i < MXN; ++i) {//预处理莫比乌斯卷积莫比乌斯
for(int j = i; j < MXN; j += i) {
mumu[j] += mu[i]*mu[j/i];
if(mumu[j] >= MOD) mumu[j] %= MOD;
}
}
for(int i = 2; i < MXN; ++i) mumu[i] = (mumu[i]+mumu[i-1]+MOD)%MOD;
}
LL solve_u(LL n) {//求莫比乌斯前缀和
if(n < MXN) return pre_mu[n];
if(mp1.count(n)) return mp1[n];
LL ans = 1;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans - (R-L+1)%MOD*solve_u(n/L)%MOD+MOD)%MOD;
}
mp1[n] = ans;
return ans;
}
LL solve_uu(LL n) {//求莫比乌斯卷积莫比乌斯前缀和
if(n < MXN) return mumu[n];
if(mp2.count(n)) return mp2[n];
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (solve_u(R)-solve_u(L-1)+MOD)%MOD*solve_u(n/L)%MOD);
if(ans >= MOD) ans %= MOD;
}
mp2[n] = (ans+MOD)%MOD;
return ans;
}
int main(int argc, char const *argv[]) {
init_rime();
scanf("%lld", &n);
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (n/L)%MOD*(n/L)%MOD*(solve_uu(R)-solve_uu(L-1)+MOD)%MOD)%MOD;
}
printf("%lld\n", (ans+MOD)%MOD);
return 0;
}
/*
void init_rime() {//因为莫比乌斯卷积莫比乌斯是积性函数,所以可以On预处理
noprime[0] = noprime[1] = 1;
mu[1] = 1;mumu[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i]=-1, mumu[i]=-2;
for(int j = 0; j < pcnt && i*pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
mu[i*pp[j]] = -mu[i];
mumu[i*pp[j]] = mumu[i]*mumu[pp[j]];
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
if((i/pp[j])%pp[j]) mumu[i*pp[j]] = mumu[i/pp[j]];
else mumu[i*pp[j]] = 0;
break;
}
}
}
for(int i = 1; i < MXN; ++i) pre_mu[i] = pre_mu[i-1] + mu[i];
for(int i = 2; i < MXN; ++i) mumu[i] = (mumu[i]+mumu[i-1]+MOD)%MOD;
}
*/
求 ∑ i = 1 n ∑ j = 1 m l c m ( i , j ) 求\sum_{i=1}^n\sum_{j=1}^{m}lcm(i,j) 求∑i=1n∑j=1mlcm(i,j)
2154 化简一下:
a
n
s
=
∑
i
=
1
n
∑
j
=
1
m
l
c
m
(
i
,
j
)
=
∑
i
=
1
n
∑
j
=
1
m
i
j
g
c
d
(
i
,
j
)
ans=\sum_{i=1}^n\sum_{j=1}^{m}lcm(i,j)=\sum_{i=1}^n\sum_{j=1}^{m}\frac{ij}{gcd(i,j)}
ans=i=1∑nj=1∑mlcm(i,j)=i=1∑nj=1∑mgcd(i,j)ij
经过一系列化简可得:
a n s = ∑ d = 1 n d ∑ k = 1 m i n ( n d , m d ) μ ( k ) × k 2 ∑ i = 1 n d k i ∑ j = 1 m k d j , O ( n ) ans=\sum_{d=1}^nd\sum_{k=1}^{min(\frac nd,\frac md)}\mu(k)\times k^2\sum_{i=1}^{\frac n{dk}}i\sum_{j=1}^{\frac m{kd}}j,O(n) ans=d=1∑ndk=1∑min(dn,dm)μ(k)×k2i=1∑dknij=1∑kdmj,O(n)
这个式子是 O ( n ) O(n) O(n),但是通过改变枚举方法可以 O ( n ) O(\sqrt n) O(n)求解。
令 T = k d , h ( n ) = ∑ i = 1 n i T=kd,h(n)=\sum_{i=1}^ni T=kd,h(n)=∑i=1ni可得:
a
n
s
=
∑
d
=
1
n
d
∑
k
=
1
n
d
μ
(
k
)
×
k
2
×
h
(
n
T
)
×
h
(
m
T
)
ans=\sum_{d=1}^nd\sum_{k=1}^{\frac nd}\mu(k)\times k^2\times h(\frac nT)\times h(\frac mT)
ans=d=1∑ndk=1∑dnμ(k)×k2×h(Tn)×h(Tm)
a
n
s
=
∑
T
=
1
n
h
(
n
T
)
h
(
m
T
)
∑
d
∣
T
d
×
μ
(
T
d
)
×
(
T
d
)
2
ans=\sum_{T=1}^nh(\frac nT)h(\frac mT)\sum_{d|T}d\times \mu(\frac Td)\times (\frac Td)^2
ans=T=1∑nh(Tn)h(Tm)d∣T∑d×μ(dT)×(dT)2
a
n
s
=
∑
T
=
1
n
h
(
n
T
)
h
(
m
T
)
∑
d
∣
T
μ
(
d
)
×
d
×
T
ans=\sum_{T=1}^nh(\frac nT)h(\frac mT)\sum_{d|T}\mu(d)\times d\times T
ans=T=1∑nh(Tn)h(Tm)d∣T∑μ(d)×d×T
然后,你发现
f
(
n
)
=
∑
d
∣
n
d
×
μ
(
d
)
f(n)=\sum_{d|n}d\times \mu(d)
f(n)=∑d∣nd×μ(d)是积性函数。为什么呢?
证明1:
i
n
v
(
x
)
=
x
−
1
inv(x)=x^{-1}
inv(x)=x−1是积性函数。
f
(
n
)
=
n
×
∑
d
∣
n
1
d
μ
(
n
d
)
=
n
×
i
n
v
∗
μ
(
n
)
f(n)=n\times \sum_{d|n}\frac 1d \mu(\frac nd)=n\times inv*\mu(n)
f(n)=n×∑d∣nd1μ(dn)=n×inv∗μ(n)
证明2:
s
(
x
)
=
x
×
μ
(
x
)
s(x)=x\times \mu(x)
s(x)=x×μ(x)是积性函数
f
(
x
)
=
s
∗
I
(
x
)
=
∑
d
∣
n
d
×
μ
(
d
)
f(x)=s*I(x)=\sum_{d|n}d\times \mu(d)
f(x)=s∗I(x)=∑d∣nd×μ(d)
很自然的可以得到:
f
(
n
)
f(n)
f(n)是一个积性函数。
当
(
a
,
b
)
=
1
(a,b)=1
(a,b)=1时,
a
a
a的因子
j
j
j和
b
b
b的因子
k
k
k当然也是互质的,然后
μ
\mu
μ是积性函数。
所以 j × μ ( j ) × k × μ ( k ) = j k × μ ( j k ) j\times \mu(j)\times k\times \mu(k)=jk\times \mu(jk) j×μ(j)×k×μ(k)=jk×μ(jk)
也就是说 f ( a ) f(a) f(a)的某一项乘上 f ( b ) f(b) f(b)的某一项一定是 f ( a b ) f(ab) f(ab)的某一项且是完备的。
f
(
1
)
=
1
,
f
(
a
)
=
1
−
a
,
a
∈
p
r
i
m
e
f(1)=1,f(a)=1-a,\;\;a\in prime
f(1)=1,f(a)=1−a,a∈prime
当
g
c
d
(
a
,
b
)
=
1
gcd(a,b)=1
gcd(a,b)=1时:
f
(
a
b
)
=
f
(
a
)
×
f
(
b
)
f(ab)=f(a)\times f(b)
f(ab)=f(a)×f(b)
当
g
c
d
(
a
,
b
)
!
=
1
gcd(a,b)!=1
gcd(a,b)!=1且
b
b
b是
a
a
a的最小素因子时:
f
(
a
b
)
=
f
(
a
)
f(ab)=f(a)
f(ab)=f(a)
为什么呢?
我们现在有:
a
a
a是
b
b
b的倍数,
f
(
a
b
)
=
∑
d
∣
(
a
b
)
d
×
μ
(
d
)
f(ab)=\sum_{d|(ab)}d\times \mu(d)
f(ab)=∑d∣(ab)d×μ(d):
如果
a
b
ab
ab的因数
d
d
d素因数分解后
b
b
b的系数小于等于
1
1
1,那么这一项一定在
f
(
a
)
f(a)
f(a)中也会出现;
如果
a
b
ab
ab的因数
d
d
d素因数分解后
b
b
b的系数大于
1
1
1,那么
μ
(
d
)
=
0
\mu(d)=0
μ(d)=0。故上式成立。
证明完 f ( n ) = ∑ d ∣ n d × μ ( d ) f(n)=\sum_{d|n}d\times \mu(d) f(n)=∑d∣nd×μ(d)是积性函数和上述式子, f ( n ) 和 f ( n ) × n f(n)和f(n)\times n f(n)和f(n)×n都可以线性筛的同时预处理出来:
void init_rime() {
noprime[0] = noprime[1] = 1;
F[1] = sum_F[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, F[i] = 1 - i;
for(int j = 0; j < pcnt && i*pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
F[i*pp[j]] = F[pp[j]]*F[i];
if(i % pp[j] == 0) {
F[i*pp[j]] = F[i];
break;
}
}
sum_F[i] = sum_F[i-1] + F[i]*i;
}
}
所以
O
(
n
)
O(\sqrt n)
O(n)即可求解:
a
n
s
=
∑
T
=
1
n
h
(
n
T
)
h
(
m
T
)
×
f
(
T
)
×
T
,
O
(
n
)
ans=\sum_{T=1}^nh(\frac nT)h(\frac mT)\times f(T)\times T,O(\sqrt n)
ans=T=1∑nh(Tn)h(Tm)×f(T)×T,O(n)
求 ∑ i = 1 n ∑ j = 1 m d ( i × j ) 求\sum_{i=1}^n\sum_{j=1}^{m}d(i\times j) 求∑i=1n∑j=1md(i×j)
3994 有一个公式:
d
(
i
×
j
)
=
∑
x
∣
i
∑
y
∣
j
[
g
c
d
(
x
,
y
)
=
1
]
d(i\times j)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]
d(i×j)=x∣i∑y∣j∑[gcd(x,y)=1]
a
n
s
=
∑
i
=
1
n
∑
x
∣
i
∑
j
=
1
m
∑
y
∣
j
[
g
c
d
(
x
,
y
)
=
1
]
ans=\sum_{i=1}^n\sum_{x|i}\sum_{j=1}^m\sum_{y|j}[gcd(x,y)=1]
ans=i=1∑nx∣i∑j=1∑my∣j∑[gcd(x,y)=1]
a
n
s
=
∑
k
=
1
n
μ
(
k
)
×
∑
x
=
1
n
k
n
k
x
×
∑
y
=
1
m
k
m
k
x
,
O
(
n
)
ans=\sum_{k=1}^n\mu(k)\times \sum_{x=1}^{\frac nk}\frac n{kx}\times \sum_{y=1}^{\frac mk}\frac m{kx},\;O(\sqrt n)
ans=k=1∑nμ(k)×x=1∑knkxn×y=1∑kmkxm,O(n)
∑
i
=
1
n
n
i
=
∑
i
=
1
n
d
(
i
)
\sum_{i=1}^n \frac ni=\sum_{i=1}^nd(i)
i=1∑nin=i=1∑nd(i)
因为
n
i
\frac ni
in表示的是
1
,
.
.
,
n
1,..,n
1,..,n有几个数是
i
i
i的倍数。当然也可以
n
n
n\sqrt n
nn预处理。
O
(
n
)
O(n)
O(n)筛法:
bool noprime[MXN];
int pp[MXN/2], pcnt, mu[MXN], pre_mu[MXN];
int a1[MXN];
LL g[MXN];
void init_prime(int MXN) {
noprime[0] = noprime[1] = mu[1] = pre_mu[1] = 1;
g[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i] = -1, g[i] = 2, a1[i] = 1;
for(int j = 0; j < pcnt && i * pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1; mu[i*pp[j]] = -mu[i];
g[i*pp[j]] = g[i] * 2;
a1[i*pp[j]] = 1;
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
a1[i*pp[j]] = a1[i] + 1;
g[i*pp[j]] = g[i]*(a1[i*pp[j]]+1)/a1[i*pp[j]];
break;
}
}
pre_mu[i] = pre_mu[i-1] + mu[i];
}
for(int i = 2; i < MXN; ++i) g[i] += g[i-1];
}
预处理
2
63
2^{63}
263范围内
d
(
n
)
d(n)
d(n)的前缀和:洛谷SP26073
小D的demo
牛客练习赛40:
x
=
∏
i
=
1
q
p
i
k
i
,
g
(
x
)
=
∑
i
=
1
q
k
i
,
g
(
1
)
=
1
x=\prod_{i=1}^qp_i^{k_i},g(x)=\sum_{i=1}^qk_i,g(1)=1
x=i=1∏qpiki,g(x)=i=1∑qki,g(1)=1
g
(
p
r
i
m
e
)
=
1
,
g
(
a
×
p
r
i
m
e
)
=
g
(
a
)
+
1
g(prime)=1,g(a\times prime)=g(a)+1
g(prime)=1,g(a×prime)=g(a)+1
所以线性筛的时候,可以顺便把
g
(
n
)
g(n)
g(n)给筛出来。
待求解式子:
∏
i
=
1
n
∏
j
=
1
m
g
(
g
c
d
(
i
,
j
)
)
\prod_{i=1}^n\prod_{j=1}^m g(gcd(i,j))
i=1∏nj=1∏mg(gcd(i,j))
提取
g
c
d
gcd
gcd:
∏
d
=
1
n
g
(
d
)
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
d
]
\prod_{d=1}^n g(d)^{\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]}
d=1∏ng(d)∑i=1n∑j=1m[gcd(i,j)=d]
反演一下指数那个式子:
∑
i
=
1
n
∑
j
=
1
m
[
g
c
d
(
i
,
j
)
=
d
]
=
∑
i
=
1
n
d
∑
j
=
1
m
d
[
g
c
d
(
i
,
j
)
=
1
]
=
∑
i
=
1
n
d
∑
j
=
1
m
d
∑
k
∣
g
c
d
(
i
,
j
)
μ
(
k
)
\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]=\sum_{i=1}^{\frac nd}\sum_{j=1}^{\frac md}[gcd(i,j)=1]=\sum_{i=1}^{\frac nd}\sum_{j=1}^{\frac md}\sum_{k|gcd(i,j)}\mu(k)
i=1∑nj=1∑m[gcd(i,j)=d]=i=1∑dnj=1∑dm[gcd(i,j)=1]=i=1∑dnj=1∑dmk∣gcd(i,j)∑μ(k)
得到这个式子:
∏
d
=
1
n
g
(
d
)
∑
k
=
1
n
d
μ
(
k
)
n
k
d
m
k
d
\prod_{d=1}^n g(d)^{\sum_{k=1}^{\frac nd}\mu(k)\frac{n}{kd}\frac{m}{kd}}
d=1∏ng(d)∑k=1dnμ(k)kdnkdm
上面这个式子可以
O
(
n
)
O(n)
O(n)求值了,但是还不够,因为这题有多组数据。
想到一个常用的优化方式:
先把指数提下来:
∏
d
=
1
n
∏
k
=
1
n
d
g
(
d
)
μ
(
k
)
n
k
d
m
k
d
\prod_{d=1}^n \prod_{k=1}^{\frac nd} g(d)^{\mu(k)\frac{n}{kd}\frac{m}{kd}}
d=1∏nk=1∏dng(d)μ(k)kdnkdm
令
T
=
k
d
T=kd
T=kd,换元求积:
∏
T
=
1
n
∏
d
∣
T
g
(
d
)
μ
(
T
d
)
n
T
m
T
\prod_{T=1}^n\prod_{d|T}g(d)^{\mu(\frac Td)\frac nT \frac mT}
T=1∏nd∣T∏g(d)μ(dT)TnTm
再变换一下:
∏
T
=
1
n
(
∏
d
∣
T
g
(
d
)
μ
(
T
d
)
)
n
T
m
T
\prod_{T=1}^n(\prod_{d|T}g(d)^{\mu(\frac Td)})^{\frac nT \frac mT}
T=1∏n(d∣T∏g(d)μ(dT))TnTm
发现这个式子
∏
d
∣
T
g
(
d
)
μ
(
T
d
)
\prod_{d|T}g(d)^{\mu(\frac Td)}
∏d∣Tg(d)μ(dT)可以
O
(
n
l
o
g
(
n
)
)
O(nlog(n))
O(nlog(n))预处理出来,
O
(
1
)
O(1)
O(1)求值。
所以最后总的复杂度是 O ( n l o g ( n ) + T n l o g ( M O D ) ) O(nlog(n)+T\sqrt n log(MOD)) O(nlog(n)+Tnlog(MOD))
代码:小D的demo
2440: [中山市选2011]完全平方数
T
=
50
T=50
T=50,求第
k
(
1
e
9
)
k(1e9)
k(1e9)个非完全平方数倍数的数字。
显然二分是一个不错的选择,当你列出求
n
n
n以内完全平方数倍数的个数的式子后,你会发现这东西不是可以用莫比乌斯函数求吗?都不需要反演的。复杂度就是枚举因数的复杂度:
O
(
n
)
O(\sqrt n)
O(n)。
HDU 6134
n
≤
1
e
6
n\leq 1e6
n≤1e6,求这个式子的值:
f
(
n
)
=
∑
d
=
1
n
μ
(
d
)
∑
x
=
1
n
d
∑
y
=
1
x
⌈
x
y
⌉
f(n)=\sum_{d=1}^n\mu(d)\sum_{x=1}^{\frac nd}\sum_{y=1}^x ⌈\frac xy⌉
f(n)=d=1∑nμ(d)x=1∑dny=1∑x⌈yx⌉
牛批网友说:
g
(
n
)
=
∑
i
=
1
n
⌈
n
i
⌉
,
ϕ
(
n
)
∑
i
=
1
n
g
(
i
)
g(n)=\sum_{i=1}^n⌈\frac ni⌉,\phi(n)\sum_{i=1}^n g(i)
g(n)=i=1∑n⌈in⌉,ϕ(n)i=1∑ng(i)
BZOJ 3529 [Sdoi2014]数表
T
(
2
e
4
)
n
,
m
(
1
e
5
)
,
A
(
1
e
9
)
T(2e4) n,m(1e5),A(1e9)
T(2e4)n,m(1e5),A(1e9),求:
∑
i
=
1
n
∑
j
=
1
m
σ
(
g
c
d
(
i
,
j
)
)
\sum_{i=1}^n\sum_{j=1}^m\sigma(gcd(i,j))
i=1∑nj=1∑mσ(gcd(i,j))
σ
(
n
)
=
∑
d
∣
n
d
;
i
f
(
σ
(
n
)
>
A
)
σ
(
n
)
=
0
;
\sigma(n)=\sum_{d|n}d;if(\sigma(n)>A)\;\sigma(n)=0;
σ(n)=d∣n∑d;if(σ(n)>A)σ(n)=0;
若不考虑
A
A
A,可以很容易得到:
∑
d
=
1
n
σ
(
d
)
∑
k
=
1
n
d
μ
(
k
)
n
k
d
n
k
d
\sum_{d=1}^n\sigma(d)\sum_{k=1}^{\frac nd}\mu(k)\frac n{kd}\frac n{kd}
d=1∑nσ(d)k=1∑dnμ(k)kdnkdn
再优化一下:
∑
T
=
1
n
n
T
m
T
∑
d
∣
T
σ
(
d
)
×
μ
(
T
d
)
\sum_{T=1}^n\frac nT\frac mT\sum_{d|T}\sigma(d)\times \mu(\frac Td)
T=1∑nTnTmd∣T∑σ(d)×μ(dT)
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MXN = 1e5 + 20;
const LL mod = 1LL<<31;
/*
T(2e4) n,m(1e5)
不考虑A的话,式子很好化简,可以sqrt(n)算答案
这题不离线下来,处理巨奇怪,也不清楚正确性
离线下来的话,主要复杂度是O(n*sqrt(n)*log(n))吧?
先把莫比乌斯函数和因数和函数给线性筛出来,离线询问按A排序
不考虑A的话,我们正常做法就是枚举gcd,用的是g(i)算答案
考虑了A,我们还要按g(i)从小到大的顺序来枚举gcd,算出卷积的贡献
因为要求卷积的前缀和,所以要用一个数据结构维护,树状数组就够了,所以复杂度多了log。
取模也可以不用取模,因为模数是1LL<<31,利用int的自然溢出就行了
*/
bool noprime[MXN];
int pp[MXN/2], pcnt;
int mu[MXN], pre_mu[MXN];
int num1[MXN];//p1^a1
int g[MXN];//因数和函数
int p[MXN];
struct lp {
int n, m, A, id;
}cw[MXN];
int ANS[MXN];
int bit[MXN];
int lowbit(int x) {return x&(-x);}
void add_bit(int x, int val, int N) {for(;x <= N; x += lowbit(x)) bit[x] = (bit[x]+val);}
int query_bit(int x) {int ans = 0; for(; x; x -= lowbit(x)) ans = (ans+bit[x]); return ans;}
void init_prime(int MXN) {
noprime[0] = noprime[1] = mu[1] = pre_mu[1] = 1;
g[1] = 1; p[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, mu[i] = -1, g[i] = 1+i, num1[i] = i;
for(int j = 0; j < pcnt && i * pp[j] < MXN; ++j) {
noprime[i*pp[j]] = 1;
mu[i*pp[j]] = -mu[i];
g[i*pp[j]] = g[i] * (1+pp[j]);
num1[i*pp[j]] = pp[j];
if(i % pp[j] == 0) {
mu[i*pp[j]] = 0;
num1[i*pp[j]] = num1[i] * pp[j];
g[i*pp[j]] = (g[i]/g[num1[i]])*(g[num1[i]]*pp[j]+1);
break;
}
}
pre_mu[i] = pre_mu[i-1] + mu[i]; p[i] = i;
}
}
bool cmp(const lp &a, const lp &b) {
return a.A < b.A;
}
bool cmp1(const int &a, const int &b) {
return g[a] < g[b];
}
int solve(int n, int m) {
int ans = 0, lst = 0, now;//除法分块的一点点优化,还是可以快一点哒
for(int L = 1, R; L <= n; L = R + 1) {
R = min(n/(n/L), m/(m/L));
now = query_bit(R);
ans = (ans+(n/L)*(m/L)*(now-lst));
lst = now;
}
return (ans);
}
int main(){
int tim; scanf("%d", &tim);
int ret = 0;
for(int i = 1; i <= tim; ++i) {
scanf("%d%d%d", &cw[i].n, &cw[i].m, &cw[i].A);
if(cw[i].n > cw[i].m) swap(cw[i].n, cw[i].m);
cw[i].id = i;
ret = max(ret, cw[i].n);
}
init_prime(ret+10);
sort(cw + 1, cw + 1 + tim, cmp);
int j = 1;
sort(p + 1, p + 1 + ret, cmp1);
for(int i = 1; i <= tim; ++i) {
for(; j <= ret && g[p[j]] <= cw[i].A; ++j) {
for(int k = p[j]; k <= ret; k += p[j]) {
if(mu[k/p[j]]) add_bit(k, g[p[j]]*mu[k/p[j]], ret);
}
}
ANS[cw[i].id] = solve(cw[i].n, cw[i].m);
}
for(int i = 1; i <= tim; ++i) {
if(ANS[i] < 0) ANS[i] += 2147483647, ++ ANS[i];
printf("%d\n", ANS[i]);
}
return 0;
}
2019西安邀请赛B. Product
显然是一个欧拉降幂套上一个莫比乌斯反演,看这数据范围还得杜教筛。
由:
∏
i
=
1
n
∏
j
=
1
n
∏
k
=
1
n
m
g
c
d
(
i
,
j
)
[
k
∣
g
c
d
(
i
,
j
)
]
\prod_{i=1}^{n}\prod_{j=1}^{n}\prod_{k=1}^{n}m^{gcd(i,j)[k|gcd(i,j)]}
i=1∏nj=1∏nk=1∏nmgcd(i,j)[k∣gcd(i,j)]
可化简为:
m
∑
i
=
1
n
∑
j
=
1
n
g
c
d
(
i
,
j
)
×
D
(
g
c
d
(
i
,
j
)
)
m^{\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)\times D(gcd(i,j))}
m∑i=1n∑j=1ngcd(i,j)×D(gcd(i,j))
意思转化为求:
∑
d
=
1
n
d
×
D
(
d
)
∑
i
=
1
n
/
d
∑
j
=
1
n
/
d
[
(
i
,
j
)
=
=
1
]
\sum_{d=1}^n d\times D(d)\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[(i,j)==1]
d=1∑nd×D(d)i=1∑n/dj=1∑n/d[(i,j)==1]
参考小清新数论那题:
∑
d
=
1
n
d
×
D
(
d
)
×
(
2
×
P
(
n
d
)
−
1
)
\sum_{d=1}^n d\times D(d)\times (2\times P(\frac nd)-1)
d=1∑nd×D(d)×(2×P(dn)−1)
又因为:
∑
i
=
1
n
i
×
D
(
i
)
=
∑
i
=
1
n
i
×
∑
d
∣
i
1
=
∑
d
=
1
n
d
×
∑
i
=
1
n
d
1
=
∑
d
×
n
d
(
n
d
+
1
)
2
\sum_{i=1}^n i\times D(i)=\sum_{i=1}^n i\times \sum_{d|i}1=\sum_{d=1}^n d\times \sum_{i=1}^\frac nd 1=\sum d\times \frac{\frac nd(\frac nd+1)}{2}
i=1∑ni×D(i)=i=1∑ni×d∣i∑1=d=1∑nd×i=1∑dn1=∑d×2dn(dn+1)
LL n, m;
int noprime[MXN], pp[MXN], pcnt;
int phi[MXN], mu[MXN], a1[MXN];
LL pre_phi[MXN], pre_g[MXN], g[MXN];
void init_prime() {//线性筛预处理
noprime[0] = noprime[1] = 1;
mu[1] = 1; phi[1] = 1; g[1] = 1;
for(int i = 2; i < MXN; ++i) {
if(!noprime[i]) pp[pcnt++] = i, phi[i] = i-1, mu[i] = -1, g[i] = 2, a1[i] = 1;
for(int j = 0; j < pcnt && pp[j] * i < MXN; ++j) {
noprime[pp[j]*i] = 1;
phi[pp[j]*i] = (pp[j]-1)*phi[i];
mu[pp[j]*i] = -mu[i];
g[i*pp[j]] = g[i] * 2;
a1[i*pp[j]] = 1;
if(i % pp[j] == 0) {
phi[pp[j]*i] = pp[j]*phi[i];
mu[pp[j]*i] = 0;
a1[i*pp[j]] = a1[i] + 1;
g[i*pp[j]] = g[i]*(a1[i*pp[j]]+1)/a1[i*pp[j]];
break;
}
}
}
}
unordered_map<LL, LL> mp1, mp2;
LL solve_g(LL n) {
if(n < MXN) return pre_g[n];
if(mp1.find(n) != mp1.end()) return mp1[n];
//if(mp1[n]) return mp1[n];
//if(mp1.count(n)) return mp1[n];
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + ((n/L)*(n/L+1)/2)%mod*((R-L+1)*(L+R)/2)%mod)%mod;
}
ans = (ans + mod) % mod;
mp1[n] = ans;
return ans;
}
LL solve_phi(LL n) {//杜教筛欧拉函数前缀和
if(n < MXN) return pre_phi[n];
if(mp2.find(n) != mp2.end()) return mp2[n];
//if(mp2[n]) return mp2[n];
//if(mp2.count(n)) return mp2[n];
LL ans = n*(n+1)/2;
for(LL L = 2, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans - (R-L+1)*solve_phi(n/L)%mod)%mod;
}
ans = (ans + mod) % mod;
mp2[n] = ans;
return ans;
}
int main() {
init_prime();
scanf("%lld%lld%lld", &n, &m, &mud);
mod = mud - 1;
for(int i = 1; i < MXN; ++i) {
pre_phi[i] = (pre_phi[i-1] + phi[i])%mod;
pre_g[i] = (pre_g[i-1] + i * g[i] % mod)%mod;
}
LL ans = 0;
for(LL L = 1, R; L <= n; L = R + 1) {
R = n/(n/L);
ans = (ans + (solve_g(R)-solve_g(L-1))%mod*(solve_phi(n/L)*2LL%mod-1LL)%mod)%mod;
}
ans = (ans+mod)%mod;
printf("%lld\n", ksm(m, ans, mud));
return 0;
}
// 当x≥ϕ(p)时,有a^x ≡ a^(x mod ϕ(p)+ϕ(p)) mod p
// gcd(a,mod)=1 a^(x) = a^(x%ϕ(mod)) % mod
南昌邀请赛和湘潭邀请赛
计蒜之道2019 1-D. 商汤AI园区的n个路口(困难)
插曲:二项式反演
F
(
n
)
=
∑
s
n
C
n
i
G
(
i
)
F(n)=\sum_{s}^{n}C_n^iG(i)
F(n)=s∑nCniG(i)
G
(
n
)
=
∑
i
=
s
n
(
−
1
)
n
−
i
C
n
i
F
(
i
)
G(n)=\sum_{i=s}^n(-1)^{n-i}C_n^iF(i)
G(n)=i=s∑n(−1)n−iCniF(i)
插曲:集合反演
m
a
x
(
S
)
=
∑
T
∈
S
(
−
1
)
∣
T
∣
−
1
×
m
i
n
(
T
)
max(S)=\sum_{T\in S}(-1)^{|T|-1}\times min(T)
max(S)=∑T∈S(−1)∣T∣−1×min(T)
m
i
n
(
S
)
=
∑
T
∈
S
(
−
1
)
∣
T
∣
−
1
×
m
a
x
(
T
)
min(S)=\sum_{T\in S}(-1)^{|T|-1}\times max(T)
min(S)=∑T∈S(−1)∣T∣−1×max(T)