莫比乌斯反演
定义
莫比乌斯反演公式:
[
n
=
1
]
=
∑
d
∣
n
μ
(
d
)
[n=1]=\underset{d|n}\sum \mu(d)
[n=1]=d∣n∑μ(d)
(可知
[
n
=
m
]
=
[
n
m
=
1
]
[n=m]=\left[\frac nm=1\right]
[n=m]=[mn=1])
其他几种莫比乌斯反演的形式:
标准形式(约数形式): f ( n ) = ∑ d ∣ n g ( d ) ⇔ g ( n ) = ∑ d ∣ n μ ( d ) f ( n d ) f(n)=\underset{d|n}\sum g(d)\Leftrightarrow g(n)=\underset{d|n}\sum\mu(d)f(\frac nd) f(n)=d∣n∑g(d)⇔g(n)=d∣n∑μ(d)f(dn),所以有一种莫反的方法是设函数。
狄利克雷卷积形式: f = g ∗ 1 ⇔ g = f ∗ μ f=g*1\Leftrightarrow g=f*\mu f=g∗1⇔g=f∗μ
倍数形式:
f
(
n
)
=
∑
n
∣
d
m
g
(
d
)
⇔
g
(
n
)
=
∑
n
∣
d
m
μ
(
d
n
)
f
(
d
)
f(n)=\overset{m}{\underset{n|d}\sum} g(d)\Leftrightarrow g(n)=\overset{m}{\underset{n|d}\sum }\mu\left(\frac dn\right)f(d)
f(n)=n∣d∑mg(d)⇔g(n)=n∣d∑mμ(nd)f(d)
现在给出莫比乌斯反演公式的证明:
对于 n = 1 n=1 n=1,显然公式成立。
对于
n
≠
1
n\neq 1
n=1,设
n
′
n'
n′为
n
n
n的所有不同质因子的连乘积。
形式化的说,有
n
n
n的标准分解式:
n
=
∏
i
=
1
s
p
i
a
i
(
p
i
−
1
<
p
i
)
n=\overset{s}{\underset{i=1}\prod}p_i^{a_i}\;\;\;(p_{i-1}<p_i)
n=i=1∏spiai(pi−1<pi),则
n
′
=
∏
i
=
1
s
p
i
(
p
i
−
1
<
p
i
)
n'=\overset{s}{\underset{i=1}\prod}p_i\;\;\;(p_{i-1}<p_i)
n′=i=1∏spi(pi−1<pi)。
显然 ∑ d ∣ n μ ( d ) = ∑ d ∣ n ′ μ ( d ) \underset{d|n}\sum\mu(d)=\underset{d|n'}\sum\mu(d) d∣n∑μ(d)=d∣n′∑μ(d),因为一旦 d d d包含相同质因子,则 μ ( d ) = 0 \mu(d)=0 μ(d)=0。
此时,假设从 n ′ n' n′中选出具有 i i i个质因子的约数 d d d,则情况数共有 ( s i ) \begin{pmatrix} s\\ i \end{pmatrix} (si)种,此时的莫比乌斯函数值是 ( − 1 ) i (-1)^i (−1)i,对答案做的贡献就是 ( − 1 ) i ( s i ) (-1)^i\begin{pmatrix} s\\ i \end{pmatrix} (−1)i(si)。枚举 d d d的质因子数目,就得到了总的贡献: ∑ i = 0 s ( − 1 ) i ( s i ) \overset{s}{\underset{i=0}\sum}(-1)^i\begin{pmatrix} s\\ i \end{pmatrix} i=0∑s(−1)i(si)
应用二项式定理:
∑
d
∣
n
μ
(
d
)
=
∑
d
∣
n
′
μ
(
d
)
=
∑
i
=
0
s
(
−
1
)
i
(
s
i
)
=
(
1
−
1
)
s
\underset{d|n}\sum\mu(d)=\underset{d|n'}\sum\mu(d)=\overset{s}{\underset{i=0}\sum}(-1)^i\begin{pmatrix} s\\ i \end{pmatrix}=(1-1)^s
d∣n∑μ(d)=d∣n′∑μ(d)=i=0∑s(−1)i(si)=(1−1)s
因此,若 n ≠ 1 n\neq 1 n=1,则 s ≥ 1 s\geq 1 s≥1,此式为 0 0 0。
QED.
对于标准形式和狄利克雷卷积形式的证明,参照狄利克雷卷积。
倍数形式的证明如下:
只证明左推右,另一方向同理:
g
(
n
)
=
∑
n
∣
d
m
μ
(
d
n
)
f
(
d
)
g(n)=\overset{m}{\underset{n|d}\sum }\mu\left(\frac dn\right)f(d)
g(n)=n∣d∑mμ(nd)f(d),将
f
f
f代入:
=
∑
n
∣
d
m
μ
(
d
n
)
∑
d
∣
d
′
m
g
(
d
′
)
=\overset{m}{\underset{n|d}\sum }\mu\left(\frac dn\right)\overset{m}{\underset{d|d'}\sum} g(d')
=n∣d∑mμ(nd)d∣d′∑mg(d′)
=
∑
n
∣
d
m
∑
d
∣
d
′
m
μ
(
d
n
)
g
(
d
′
)
=\overset{m}{\underset{n|d}\sum }\overset{m}{\underset{d|d'}\sum}\mu\left(\frac dn\right) g(d')
=n∣d∑md∣d′∑mμ(nd)g(d′)
n ∣ d , d ∣ d ′ n|d,d|d' n∣d,d∣d′,相当于先枚举 n ∣ d ′ n|d' n∣d′,再枚举 n ∣ d ∣ d ′ n|d|d' n∣d∣d′:
=
∑
n
∣
d
′
m
∑
n
∣
d
∣
d
′
m
μ
(
d
n
)
g
(
d
′
)
=\overset{m}{\underset{n|d'}\sum }\overset{m}{\underset{n|d|d'}\sum}\mu\left(\frac dn\right) g(d')
=n∣d′∑mn∣d∣d′∑mμ(nd)g(d′)
=
∑
n
∣
d
′
m
∑
d
∣
d
′
n
m
/
n
μ
(
d
)
g
(
d
′
)
=\overset{m}{\underset{n|d'}\sum }\overset{m/n}{\underset{d|\frac{d'}n}\sum}\mu\left(d\right) g(d')
=n∣d′∑md∣nd′∑m/nμ(d)g(d′)
=
∑
n
∣
d
′
m
g
(
d
′
)
∑
d
∣
d
′
n
m
/
n
μ
(
d
)
=\overset{m}{\underset{n|d'}\sum } g(d')\overset{m/n}{\underset{d|\frac{d'}n}\sum}\mu\left(d\right)
=n∣d′∑mg(d′)d∣nd′∑m/nμ(d)
因为 m / n ≥ d ′ n m/n\geq\frac {d'}n m/n≥nd′:
=
∑
n
∣
d
′
m
g
(
d
′
)
∑
d
∣
d
′
n
μ
(
d
)
=\overset{m}{\underset{n|d'}\sum } g(d')\overset{}{\underset{d|\frac{d'}n}\sum}\mu\left(d\right)
=n∣d′∑mg(d′)d∣nd′∑μ(d)
=
∑
n
∣
d
′
m
g
(
d
′
)
[
d
′
n
=
1
]
=\overset{m}{\underset{n|d'}\sum } g(d')\left[\frac {d'}n=1\right]
=n∣d′∑mg(d′)[nd′=1]
=
∑
n
∣
d
′
m
g
(
d
′
)
[
d
′
=
n
]
=\overset{m}{\underset{n|d'}\sum } g(d')\left[ {d'}=n\right]
=n∣d′∑mg(d′)[d′=n]
=
g
(
n
)
=g(n)
=g(n)
QED.
简单莫比乌斯反演
题目要求:
∑
i
=
1
a
∑
j
=
1
b
[
gcd
(
i
,
j
)
=
d
]
\overset{a}{\underset{i=1}\sum}\overset{b}{\underset{j=1}\sum}[\gcd(i,j)=d]
i=1∑aj=1∑b[gcd(i,j)=d]
由于个人习惯,替换一下符号,设
n
≤
m
n\leq m
n≤m。
∑
i
=
1
n
∑
j
=
1
m
[
gcd
(
i
,
j
)
=
g
]
\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}[\gcd(i,j)=g]
i=1∑nj=1∑m[gcd(i,j)=g]
做和式的变换:
替换指标变量,同时改变枚举范围:
∑
i
=
1
n
∑
j
=
1
m
[
gcd
(
i
,
j
)
=
g
]
=
∑
i
=
1
⌊
n
g
⌋
∑
j
=
1
⌊
m
g
⌋
[
gcd
(
i
,
j
)
=
1
]
\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}[\gcd(i,j)=g]=\overset{\left\lfloor\frac ng\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac mg\right\rfloor}{\underset{j=1}\sum}[\gcd(i,j)=1]
i=1∑nj=1∑m[gcd(i,j)=g]=i=1∑⌊gn⌋j=1∑⌊gm⌋[gcd(i,j)=1]
带入莫比乌斯反演公式:
=
∑
i
=
1
⌊
n
g
⌋
∑
j
=
1
⌊
m
g
⌋
∑
d
∣
gcd
(
i
,
j
)
μ
(
d
)
=\overset{\left\lfloor\frac ng\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac mg\right\rfloor}{\underset{j=1}\sum}\underset{d|\gcd(i,j)}\sum\mu(d)
=i=1∑⌊gn⌋j=1∑⌊gm⌋d∣gcd(i,j)∑μ(d)
替换条件式:
=
∑
i
=
1
⌊
n
g
⌋
∑
j
=
1
⌊
m
g
⌋
∑
d
=
1
n
[
d
∣
i
]
[
d
∣
j
]
μ
(
d
)
=\overset{\left\lfloor\frac ng\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac mg\right\rfloor}{\underset{j=1}\sum}\overset{n}{\underset{d=1}\sum}[d|i][d|j]\mu(d)
=i=1∑⌊gn⌋j=1∑⌊gm⌋d=1∑n[d∣i][d∣j]μ(d)
分离变量:
=
∑
d
=
1
n
μ
(
d
)
∑
i
=
1
⌊
n
g
⌋
[
d
∣
i
]
∑
j
=
1
⌊
m
g
⌋
[
d
∣
j
]
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{\left\lfloor\frac ng\right\rfloor}{\underset{i=1}\sum}[d|i]\overset{\left\lfloor\frac mg\right\rfloor}{\underset{j=1}\sum}[d|j]
=d=1∑nμ(d)i=1∑⌊gn⌋[d∣i]j=1∑⌊gm⌋[d∣j]
注意到有: ∑ i = 1 n [ d ∣ i ] = ⌊ n d ⌋ \overset{n}{\underset{i=1}\sum}[d|i]=\left\lfloor\frac n d\right\rfloor i=1∑n[d∣i]=⌊dn⌋,因此:
= ∑ d = 1 n μ ( d ) ⌊ n d g ⌋ ⌊ m d g ⌋ =\overset{n}{\underset{d=1}\sum}\mu(d)\left\lfloor\frac n {dg}\right\rfloor\left\lfloor\frac m {dg}\right\rfloor =d=1∑nμ(d)⌊dgn⌋⌊dgm⌋
于是可以整除分块。
#include<iostream>
#include<vector>
using namespace std;
bool vis[50005];
vector<int> sta;
int mu[50005],pre[50005];
void Euler() {
mu[1]=1;
for(int i=2;i<=50000;i++) {
if(!vis[i])
sta.push_back(i),mu[i]=-1;
for(auto&j:sta) {
if(i*j>50000) break;
int m=i*j;
vis[m]=true;
if(i%j) mu[m]=-mu[i];
else break;
}
}
}
inline int sum(int l,int r) {
return pre[r]-pre[l-1];
}
int calc(int n,int m) {
int ans=0;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),
ans+=(n/l)*(m/l)*sum(l,r);
return ans;
}
int main() {
Euler();
for(int i=1;i<=50000;i++) pre[i]=pre[i-1]+mu[i];
// for(int i=1;i<=10;i++)
// cout<<mu[i]<<' ';
int T;
cin>>T;
while(T--) {
int a,b,d;
cin>>a>>b>>d;
if(a>b) swap(a,b);
cout<<calc(a/d,b/d)<<endl;
}
}
莫比乌斯反演技巧
假设n ≤ m ≤ K \leq m\leq K ≤m≤K
枚举倍数与换元
题目要求:
∑
p
∈
P
∑
i
=
1
N
∑
j
=
1
M
[
gcd
(
i
,
j
)
=
p
]
\underset{p\in\mathbb{P}}\sum\overset{N}{\underset{i=1}\sum}\overset{M}{\underset{j=1}\sum}[\gcd(i,j)=p]
p∈P∑i=1∑Nj=1∑M[gcd(i,j)=p]
改一下字母:
∑
p
∈
P
∑
i
=
1
n
∑
j
=
1
m
[
gcd
(
i
,
j
)
=
p
]
\underset{p\in\mathbb{P}}\sum\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}[\gcd(i,j)=p]
p∈P∑i=1∑nj=1∑m[gcd(i,j)=p]
替换指标变量:
=
∑
p
∈
P
∑
i
=
1
⌊
n
p
⌋
∑
j
=
1
⌊
m
p
⌋
[
gcd
(
i
,
j
)
=
1
]
=\underset{p\in\mathbb{P}}\sum\overset{\left\lfloor\frac n p\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac m p\right\rfloor}{\underset{j=1}\sum}[\gcd(i,j)=1]
=p∈P∑i=1∑⌊pn⌋j=1∑⌊pm⌋[gcd(i,j)=1]
枚举上界和变量有关,一看就不是什么好事。…但是暂时替换不了。
带入反演公式:
=
∑
p
∈
P
∑
i
=
1
⌊
n
p
⌋
∑
j
=
1
⌊
m
p
⌋
∑
d
∣
gcd
(
i
,
j
)
μ
(
d
)
=\underset{p\in\mathbb{P}}\sum\overset{\left\lfloor\frac n p\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac m p\right\rfloor}{\underset{j=1}\sum}\underset{d|\gcd(i,j)}\sum\mu(d)
=p∈P∑i=1∑⌊pn⌋j=1∑⌊pm⌋d∣gcd(i,j)∑μ(d)
替换条件式:
=
∑
p
∈
P
∑
i
=
1
⌊
n
p
⌋
∑
j
=
1
⌊
m
p
⌋
∑
d
=
1
⌊
n
p
⌋
[
d
∣
i
]
[
d
∣
j
]
μ
(
d
)
=\underset{p\in\mathbb{P}}\sum\overset{\left\lfloor\frac n p\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac m p\right\rfloor}{\underset{j=1}\sum}\overset{\left\lfloor\frac n p\right\rfloor}{\underset{d=1}\sum}[d|i][d|j]\mu(d)
=p∈P∑i=1∑⌊pn⌋j=1∑⌊pm⌋d=1∑⌊pn⌋[d∣i][d∣j]μ(d)
事实上枚举范围还可以扩大:
=
∑
p
∈
P
∑
i
=
1
⌊
n
p
⌋
∑
j
=
1
⌊
m
p
⌋
∑
d
=
1
n
[
d
∣
i
]
[
d
∣
j
]
μ
(
d
)
=\underset{p\in\mathbb{P}}\sum\overset{\left\lfloor\frac n p\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac m p\right\rfloor}{\underset{j=1}\sum}\overset{n}{\underset{d=1}\sum}[d|i][d|j]\mu(d)
=p∈P∑i=1∑⌊pn⌋j=1∑⌊pm⌋d=1∑n[d∣i][d∣j]μ(d)
这样最后一个和式的枚举范围就和变量无关了。
改变求和次序,分离变量:
=
∑
d
=
1
n
μ
(
d
)
∑
p
∈
P
∑
i
=
1
⌊
n
p
⌋
[
d
∣
i
]
∑
j
=
1
⌊
m
p
⌋
[
d
∣
j
]
=\overset{n}{\underset{d=1}\sum}\mu(d)\underset{p\in\mathbb{P}}\sum\overset{\left\lfloor\frac n p\right\rfloor}{\underset{i=1}\sum}[d|i]\overset{\left\lfloor\frac m p\right\rfloor}{\underset{j=1}\sum}[d|j]
=d=1∑nμ(d)p∈P∑i=1∑⌊pn⌋[d∣i]j=1∑⌊pm⌋[d∣j]
则有:
=
∑
d
=
1
n
μ
(
d
)
∑
p
∈
P
⌊
n
d
⋅
p
⌋
⌊
m
d
⋅
p
⌋
=\overset{n}{\underset{d=1}\sum}\mu(d)\underset{p\in\mathbb{P}}\sum{\left\lfloor\frac n{d\cdot p}\right\rfloor}\left\lfloor\frac m {d\cdot p}\right\rfloor
=d=1∑nμ(d)p∈P∑⌊d⋅pn⌋⌊d⋅pm⌋
此时看似是无法处理了,但是敏锐的我们发现了问题:
式子中限制了
d
⋅
p
≤
n
d\cdot p\leq n
d⋅p≤n,则我们注意到,在
[
n
]
[n]
[n]中,乘积
≤
n
\leq n
≤n的二元组数量应该为
∑
n
i
=
1
d
(
i
)
=
O
(
n
ln
n
)
\underset{i=1}{\overset n\sum}d(i)=O(n\ln n)
i=1∑nd(i)=O(nlnn)。
也就是说,对于每一次询问枚举所有的
d
,
p
d,p
d,p是不值得的。应当有一种方法预处理所有的
d
,
p
d,p
d,p对答案的贡献。
具体来说,如果 d 1 p 1 = d 2 p 2 d_1p_1=d_2p_2 d1p1=d2p2,那么他们对向下取整部分的贡献是相同的。这一部分是与 n , m n,m n,m有关的。还有一部分是他们各自的贡献,例如 μ ( d ) , [ p ∈ P ] \mu(d),[p\in\mathbb{P}] μ(d),[p∈P],这一部分我们可以预处理出来乘到前者上面,则只需要计算第一部分的贡献即可。
具体来说:
-
换元,令 T = d ⋅ p T=d\cdot p T=d⋅p,注意到 ⌊ n T ⌋ \left\lfloor\frac nT\right\rfloor ⌊Tn⌋,因此 T ≤ n T\leq n T≤n:
= ∑ T = 1 n ⌊ n T ⌋ ⌊ m T ⌋ ∑ d = 1 n [ d ∣ T ] μ ( d ) [ T d ∈ P ] =\overset{n}{\underset{T=1}\sum}{\left\lfloor\frac n{T}\right\rfloor}\left\lfloor\frac m {T}\right\rfloor\overset{n}{\underset{d=1}\sum}[d|T]\mu(d)[\frac Td\in\mathbb{P}] =T=1∑n⌊Tn⌋⌊Tm⌋d=1∑n[d∣T]μ(d)[dT∈P]
似乎无法可想? -
换元,令 T = d ⋅ p T=d\cdot p T=d⋅p,注意到 ⌊ n T ⌋ \left\lfloor\frac nT\right\rfloor ⌊Tn⌋,因此 T ≤ n T\leq n T≤n:
=
∑
T
=
1
n
⌊
n
T
⌋
⌊
m
T
⌋
∑
p
∈
P
[
p
∣
T
]
μ
(
T
p
)
=\overset{n}{\underset{T=1}\sum}\left\lfloor\frac n{T}\right\rfloor\left\lfloor\frac m {T}\right\rfloor\overset{}{\underset{p\in\mathbb{P}}\sum}[p|T]\mu(\frac T p)
=T=1∑n⌊Tn⌋⌊Tm⌋p∈P∑[p∣T]μ(pT)
这个形式看起来很是很不可做?事实上我们仔细观察右边,已经可做了。
设函数 F ( x ) = ∑ p ∈ P [ p ∣ x ] μ ( x p ) F(x)=\underset{p\in\mathbb{P}}\sum[p|x]\mu\left(\frac x p\right) F(x)=p∈P∑[p∣x]μ(px),我们注意到, p p p是固定的一些数,而 p p p只对它的倍数做贡献,因而我们枚举 p p p,枚举 p p p的倍数,要枚举的次数是 O ( ∑ p ∈ P n p ) O\left(\underset{p\in\mathbb{P}}\sum\frac n p\right) O(p∈P∑pn),根据埃氏筛时间复杂度的证明, O ( ∑ p ∈ P n p ) = O ( n ln ln n ) O\left(\underset{p\in\mathbb{P}}\sum\frac n p\right)=O(n\ln\ln n) O(p∈P∑pn)=O(nlnlnn)。
可以通过该题。
#include<iostream>
#include<vector>
using namespace std;
const int N=1e7;
bool vis[N+5];
vector<int> sta;
int mu[N+5];
long long f[N+5],pre[N+5];
void Euler() {
mu[1]=1;
for(int i=2;i<=N;i++) {
if(!vis[i])
sta.push_back(i),mu[i]=-1;
for(auto&j:sta) {
if(i*j>N) break;
int m=i*j;
vis[m]=true;
if(i%j) mu[m]=-mu[i];
else break;
}
}
for(auto&j:sta)
for(int i=j;i<=N;i+=j)
f[i]+=mu[i/j];
for(int i=1;i<=N;i++)
pre[i]=pre[i-1]+f[i];
}
long long calc(int n,int m) {
long long ans=0;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),
ans+=(pre[r]-pre[l-1])*(n/l)*(m/l);
return ans;
}
int main() {
Euler();
int T;
cin>>T;
while(T--) {
int n,m;
cin>>n>>m;
if(n>m) swap(n,m);
cout<<calc(n,m)<<endl;
}
}
其实第一种换元也是可以做的,只需要枚举 d d d的质数倍就好了。其实换元只是为了看起来清晰一点。
容斥
大矩形表示
(
b
,
d
)
(b,d)
(b,d)范围内的答案,再减去绿色、橙色范围内的答案,最后加上紫色的答案,就得到了红色阴影部分的答案,也就是要求的答案。因此做一个二维容斥。
#include<iostream>
#include<vector>
using namespace std;
bool vis[50005];
vector<int> sta;
int mu[50005],pre[50005];
void Euler() {
mu[1]=1;
for(int i=2;i<=50000;i++) {
if(!vis[i])
sta.push_back(i),mu[i]=-1;
for(auto&j:sta) {
if(i*j>50000) break;
int m=i*j;
vis[m]=true;
if(i%j) mu[m]=-mu[i];
else break;
}
}
}
inline int sum(int l,int r) {
return pre[r]-pre[l-1];
}
int calc(int n,int m) {
if(n>m) swap(n,m);
int ans=0;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),
ans+=(n/l)*(m/l)*sum(l,r);
return ans;
}
int main() {
Euler();
for(int i=1;i<=50000;i++) pre[i]=pre[i-1]+mu[i];
// for(int i=1;i<=10;i++)
// cout<<mu[i]<<' ';
int T;
cin>>T;
while(T--) {
int a,b,c,d,k;
cin>>a>>b>>c>>d>>k;
cout<<calc(b/k,d/k)-calc((a-1)/k,d/k)-calc(b/k,(c-1)/k)+calc((a-1)/k,(c-1)/k)<<endl;
}
}
约数个数函数(一)
题目要求:
∑
i
=
1
n
∑
j
=
1
m
d
(
i
⋅
j
)
\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}d(i\cdot j)
i=1∑nj=1∑md(i⋅j)
其中 d d d表述约数个数函数。
事实上必须知道约数个数函数有性质:
d
(
i
⋅
j
)
=
∑
x
∣
i
∑
y
∣
j
[
x
⊥
y
]
d(i\cdot j)={\underset{x|i}\sum}\overset{}{\underset{y|j}\sum}[x\perp y]
d(i⋅j)=x∣i∑y∣j∑[x⊥y]
这里贴出证明:
让大家喜闻乐见的证明(推荐)
接近本质的证明(不完善,不推荐)
因此可以反演:
∑
i
=
1
n
∑
j
=
1
m
d
(
i
⋅
j
)
\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}d(i\cdot j)
i=1∑nj=1∑md(i⋅j)
=
∑
i
=
1
n
∑
j
=
1
m
∑
x
∣
i
∑
y
∣
j
[
x
⊥
y
]
=\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}{\underset{x|i}\sum}\overset{}{\underset{y|j}\sum}[x\perp y]
=i=1∑nj=1∑mx∣i∑y∣j∑[x⊥y]
=
∑
i
=
1
n
∑
j
=
1
m
∑
x
=
1
n
[
x
∣
i
]
∑
y
=
1
m
[
y
∣
j
]
[
x
⊥
y
]
=\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}\overset{n}{\underset{x=1}\sum}[x|i]\overset{m}{\underset{y=1}\sum}[y|j][x\perp y]
=i=1∑nj=1∑mx=1∑n[x∣i]y=1∑m[y∣j][x⊥y]
=
∑
i
=
1
n
∑
j
=
1
m
∑
x
=
1
n
[
x
∣
i
]
∑
y
=
1
m
[
y
∣
j
]
∑
d
=
1
n
[
d
∣
x
]
[
d
∣
y
]
μ
(
d
)
=\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}\overset{n}{\underset{x=1}\sum}[x|i]\overset{m}{\underset{y=1}\sum}[y|j]\overset{n}{\underset{d=1}\sum}[d|x][d|y]\mu(d)
=i=1∑nj=1∑mx=1∑n[x∣i]y=1∑m[y∣j]d=1∑n[d∣x][d∣y]μ(d)
=
∑
d
=
1
n
μ
(
d
)
∑
x
=
1
n
[
d
∣
x
]
∑
y
=
1
m
[
d
∣
y
]
∑
i
=
1
n
[
x
∣
i
]
∑
j
=
1
m
[
y
∣
j
]
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{n}{\underset{x=1}\sum}[d|x]\overset{m}{\underset{y=1}\sum}[d|y]\overset{n}{\underset{i=1}\sum}[x|i]\overset{m}{\underset{j=1}\sum}[y|j]
=d=1∑nμ(d)x=1∑n[d∣x]y=1∑m[d∣y]i=1∑n[x∣i]j=1∑m[y∣j]
=
∑
d
=
1
n
μ
(
d
)
∑
x
=
1
n
[
d
∣
x
]
⌊
n
x
⌋
∑
y
=
1
m
[
d
∣
y
]
⌊
n
y
⌋
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{n}{\underset{x=1}\sum}[d|x]\left\lfloor\frac nx\right\rfloor\overset{m}{\underset{y=1}\sum}[d|y]\left\lfloor\frac ny\right\rfloor
=d=1∑nμ(d)x=1∑n[d∣x]⌊xn⌋y=1∑m[d∣y]⌊yn⌋
注意到此时已经无法直接把
x
,
y
x,y
x,y化掉了,因为
x
,
y
x,y
x,y不仅缀有艾弗森括号,还有其他式子。此时我们直接令
x
=
d
⋅
x
,
y
=
d
⋅
y
x=d\cdot x,y=d\cdot y
x=d⋅x,y=d⋅y,另外注意到此时我们的枚举上界是可以不变的,因为一旦
d
x
>
n
dx>n
dx>n,则和式的值为
0
0
0,因此没必要做限制(虽然本题不影响):
=
∑
d
=
1
n
μ
(
d
)
∑
x
=
1
n
⌊
n
d
x
⌋
∑
y
=
1
m
⌊
m
d
y
⌋
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{n}{\underset{x=1}\sum}\left\lfloor\frac n{dx}\right\rfloor\overset{m}{\underset{y=1}\sum}\left\lfloor\frac m{dy}\right\rfloor
=d=1∑nμ(d)x=1∑n⌊dxn⌋y=1∑m⌊dym⌋
发现已经可做了,因为整除定理表明
⌊
a
b
c
⌋
=
⌊
⌊
a
b
⌋
c
⌋
\left\lfloor\frac a {bc}\right\rfloor=\left\lfloor\frac{\lfloor\frac a b\rfloor}c\right\rfloor
⌊bca⌋=⌊c⌊ba⌋⌋:
设
F
(
x
)
=
∑
i
=
1
x
⌊
x
i
⌋
F(x)=\overset{x}{\underset{i=1}\sum}\left\lfloor\frac x{i}\right\rfloor
F(x)=i=1∑x⌊ix⌋,
F
(
x
)
F(x)
F(x)可以直接
O
(
n
n
)
O(n\sqrt n)
O(nn)预处理出来,这已经可以做了。当然事实上,接下来你会看到,这个东西可以
O
(
n
)
O(n)
O(n)预处理。
因此最后:
=
∑
d
=
1
n
μ
(
d
)
F
(
⌊
n
d
⌋
)
F
(
⌊
m
d
⌋
)
=\overset{n}{\underset{d=1}\sum}\mu(d)F\left(\left\lfloor\frac n{d}\right\rfloor\right)F\left(\left\lfloor\frac m{d}\right\rfloor\right)
=d=1∑nμ(d)F(⌊dn⌋)F(⌊dm⌋)
#include<iostream>
#include<vector>
using namespace std;
const int N=5e4;
bool vis[N+5];
vector<int> sta;
int mu[N+5];
long long pre[N+5];
long long F[N+5];
//F[i]表示∑x:[i/x]
void Euler() {
mu[1]=1;
for(int i=2;i<=N;i++) {
if(!vis[i])
sta.push_back(i),mu[i]=-1;
for(auto&j:sta) {
if(i*j>N) break;
int m=i*j;
vis[m]=true;
if(i%j) mu[m]=-mu[i];
else break;
}
}
for(int i=1;i<=N;i++) pre[i]=pre[i-1]+mu[i];
for(int i=1;i<=N;i++)
for(int l=1,r;l<=i;l=r+1)
r=(i/(i/l)),
F[i]+=(long long)(i/l)*(r-l+1);
}
long long calc(int n,int m) {
long long ans=0;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),
ans+=(pre[r]-pre[l-1])*F[n/l]*F[m/l];
return ans;
}
int main() {
Euler();
int T;
cin>>T;
while(T--) {
int n,m;
cin>>n>>m;
if(n>m) swap(n,m);
cout<<calc(n,m)<<endl;
}
}
枚举gcd,分块套分块,小学奥数
题目要求:
∑
i
=
1
n
∑
j
=
1
m
l
c
m
(
i
,
j
)
\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}lcm(i,j)
i=1∑nj=1∑mlcm(i,j)
首先肯定把
l
c
m
lcm
lcm搞成
g
c
d
gcd
gcd:
=
∑
i
=
1
n
∑
j
=
1
m
i
⋅
j
gcd
(
i
,
j
)
=\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}\frac {i\cdot j}{\gcd(i,j)}
=i=1∑nj=1∑mgcd(i,j)i⋅j
这个式子似乎已经很难做了,又没法换元,因此我们考虑直接枚举
g
c
d
gcd
gcd,暴力搞一下这个式子:
=
∑
g
=
1
n
∑
i
=
1
n
∑
j
=
1
m
i
⋅
j
g
[
gcd
(
i
,
j
)
=
g
]
=\overset{n}{\underset{g=1}\sum}\overset{n}{\underset{i=1}\sum}\overset{m}{\underset{j=1}\sum}\frac {i\cdot j}{g}[\gcd(i,j)=g]
=g=1∑ni=1∑nj=1∑mgi⋅j[gcd(i,j)=g]
=
∑
g
=
1
n
∑
i
=
1
⌊
n
g
⌋
∑
j
=
1
⌊
m
g
⌋
i
⋅
j
⋅
g
[
gcd
(
i
,
j
)
=
1
]
=\overset{n}{\underset{g=1}\sum}\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac m{g}\right\rfloor}{\underset{j=1}\sum}i\cdot j\cdot g[\gcd(i,j)=1]
=g=1∑ni=1∑⌊gn⌋j=1∑⌊gm⌋i⋅j⋅g[gcd(i,j)=1]
=
∑
g
=
1
n
g
∑
i
=
1
⌊
n
g
⌋
i
∑
j
=
1
⌊
m
g
⌋
j
[
gcd
(
i
,
j
)
=
1
]
=\overset{n}{\underset{g=1}\sum}g\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}i\overset{\left\lfloor\frac m{g}\right\rfloor}{\underset{j=1}\sum}j[\gcd(i,j)=1]
=g=1∑ngi=1∑⌊gn⌋ij=1∑⌊gm⌋j[gcd(i,j)=1]
=
∑
g
=
1
n
g
∑
i
=
1
⌊
n
g
⌋
i
∑
j
=
1
⌊
m
g
⌋
j
∑
d
=
1
⌊
n
g
⌋
[
d
∣
i
]
[
d
∣
j
]
μ
(
d
)
=\overset{n}{\underset{g=1}\sum}g\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}i\overset{\left\lfloor\frac m{g}\right\rfloor}{\underset{j=1}\sum}j\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{d=1}\sum}[d|i][d|j]\mu(d)
=g=1∑ngi=1∑⌊gn⌋ij=1∑⌊gm⌋jd=1∑⌊gn⌋[d∣i][d∣j]μ(d)
=
∑
g
=
1
n
g
∑
i
=
1
⌊
n
g
⌋
i
∑
j
=
1
⌊
m
g
⌋
j
∑
d
=
1
n
[
d
∣
i
]
[
d
∣
j
]
μ
(
d
)
=\overset{n}{\underset{g=1}\sum}g\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}i\overset{\left\lfloor\frac m{g}\right\rfloor}{\underset{j=1}\sum}j\overset{n}{\underset{d=1}\sum}[d|i][d|j]\mu(d)
=g=1∑ngi=1∑⌊gn⌋ij=1∑⌊gm⌋jd=1∑n[d∣i][d∣j]μ(d)
=
∑
d
=
1
n
μ
(
d
)
∑
g
=
1
n
g
∑
i
=
1
⌊
n
g
⌋
[
d
∣
i
]
i
∑
j
=
1
⌊
m
g
⌋
[
d
∣
j
]
j
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{n}{\underset{g=1}\sum}g\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}[d|i]i\overset{\left\lfloor\frac m{g}\right\rfloor}{\underset{j=1}\sum}[d|j]j
=d=1∑nμ(d)g=1∑ngi=1∑⌊gn⌋[d∣i]ij=1∑⌊gm⌋[d∣j]j
=
∑
d
=
1
n
μ
(
d
)
∑
g
=
1
n
g
∑
i
=
1
⌊
n
d
g
⌋
i
⋅
d
∑
j
=
1
⌊
m
d
g
⌋
j
⋅
d
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{n}{\underset{g=1}\sum}g\overset{\left\lfloor\frac n{dg}\right\rfloor}{\underset{i=1}\sum}i\cdot d\overset{\left\lfloor\frac m{dg}\right\rfloor}{\underset{j=1}\sum}j\cdot d
=d=1∑nμ(d)g=1∑ngi=1∑⌊dgn⌋i⋅dj=1∑⌊dgm⌋j⋅d
=
∑
d
=
1
n
μ
(
d
)
d
2
∑
g
=
1
n
g
∑
i
=
1
⌊
n
d
g
⌋
i
∑
j
=
1
⌊
m
d
g
⌋
j
=\overset{n}{\underset{d=1}\sum}\mu(d)d^2\overset{n}{\underset{g=1}\sum}g\overset{\left\lfloor\frac n{dg}\right\rfloor}{\underset{i=1}\sum}i\overset{\left\lfloor\frac m{dg}\right\rfloor}{\underset{j=1}\sum}j
=d=1∑nμ(d)d2g=1∑ngi=1∑⌊dgn⌋ij=1∑⌊dgm⌋j
这时候已经可以做了,不是吗?
外层对
d
d
d分块,内层对
g
g
g分块,分块套分块,时间复杂度
O
(
n
)
O(n)
O(n)。
最后两个和式是等差数列,平方数的和,小学奥数是可以速算的,再不济可以打表,对不对?
#include<iostream>
#include<vector>
using namespace std;
const int N=1e7;
bool vis[N+5];
vector<int> sta;
int mu[N+5];
void Euler() {
mu[1]=1;
for(int i=2;i<=N;i++) {
if(!vis[i])
sta.push_back(i),mu[i]=-1;
for(auto&j:sta) {
if(i*j>N) break;
int m=i*j;
vis[m]=true;
if(i%j) mu[m]=-mu[i];
else break;
}
}
}
const long long mod=20101009;
long long f[N+5],F[N+5];
long long sum(int n,int m) {
long long ans=0;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),
(ans+=((f[r]-f[l-1])%mod+mod)*f[n/l]%mod*f[m/l])%=mod;
return ans;
}
long long calc(int n,int m) {
long long ans=0;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),
(ans+=((F[r]-F[l-1])%mod+mod)*sum(n/l,m/l))%=mod;
return ans;
}
int main() {
int n,m;
Euler();
for(int i=1;i<=N;i++)
f[i]=((long long)i*(i+1)>>1)%mod;
for(int i=1;i<=N;i++)
F[i]=(((long long)i*i*mu[i]+F[i-1])%mod+mod)%mod;
cin>>n>>m;
if(n>m) swap(n,m);
cout<<calc(n,m);
}
后记:
当然还有一种做法是先对 g g g分块,再对 d d d分块,也是可以做的,不过我们能够让 d d d成为分块的外层,还归功于我们扩大了枚举范围
和式与prod式的变换
题目要求:
∏
i
=
1
n
∏
j
=
1
m
f
(
g
c
d
(
i
,
j
)
)
\overset{n}{\underset{i=1}\prod}\overset{m}{\underset{j=1}\prod}f\left(gcd(i,j)\right)
i=1∏nj=1∏mf(gcd(i,j))
prod式应该怎么变换的?其实与和式的变换差不多:
- prod式满足交换律。
- 满足结合律: ∏ i f ( i ) ∏ j g ( j ) ∏ k h ( k ) = ∏ i f ( i ) ( ∏ j g ( j ) ∏ k h ( k ) ) \underset{i}\prod f(i)\underset{j}\prod g(j)\underset{k}\prod h(k)=\underset{i}\prod f(i)\left(\underset{j}\prod g(j)\underset{k}\prod h(k)\right) i∏f(i)j∏g(j)k∏h(k)=i∏f(i)(j∏g(j)k∏h(k))
- 满足分配律:
∏
i
=
1
n
c
a
i
=
c
n
∏
i
=
1
n
a
i
\overset{n}{\underset{i=1}\prod}ca_i=c^n\overset{n}{\underset{i=1}\prod}a_i
i=1∏ncai=cni=1∏nai
至于加法…好像没法分配 - 指数: ( ∏ i f ( i ) ∏ j g ( j ) ) k = ( ∏ i f ( i ) ) k ( ∏ j g ( j ) ) k = ∏ i f k ( i ) ∏ j g k ( j ) \left(\underset{i}\prod f(i)\underset{j}\prod g(j)\right)^k=\left(\underset{i}\prod f(i)\right)^k\left(\underset{j}\prod g(j)\right)^k=\underset{i}\prod f^k(i)\underset{j}\prod g^k(j) (i∏f(i)j∏g(j))k=(i∏f(i))k(j∏g(j))k=i∏fk(i)j∏gk(j)
- 任何大型运算符都可以替换条件式,只不过细节不同,求和符号条件式转艾弗森写在系数上: ∑ d ∣ n f ( d ) = ∑ [ d ∣ n ] f ( d ) \underset{d|n}\sum f(d)=\sum[d|n]f(d) d∣n∑f(d)=∑[d∣n]f(d),而连乘积符号条件式转艾弗森写在指数上: ∏ d ∣ n f ( d ) = ∏ f ( d ) [ d ∣ n ] \underset{d|n}\prod f(d)=\prod f(d)^{[d|n]} d∣n∏f(d)=∏f(d)[d∣n]
- 任何大型运算符都可以替换指标变量
- 改变运算次序: ∏ i ∏ j f ( i ) g ( j ) = ∏ j ∏ i f ( i ) g ( j ) \underset{i}\prod\underset{j}\prod f(i)g(j)=\underset{j}\prod\underset{i}\prod f(i)g(j) i∏j∏f(i)g(j)=j∏i∏f(i)g(j)
- 分离变量: ∏ i ∏ j f ( i ) g ( j ) = ∏ i f ( i ) ∏ j g ( j ) \underset{i}\prod\underset{j}\prod f(i)g(j)=\underset{i}\prod f(i)\underset{j}\prod g(j) i∏j∏f(i)g(j)=i∏f(i)j∏g(j)
- 求和符号与连乘积符号:
∏
i
∏
j
f
(
i
)
g
(
j
)
=
∏
i
f
(
i
)
∑
j
g
(
j
)
\underset{i}\prod\underset{j}\prod f(i)^{g(j)}=\underset{i}\prod f(i)^{\underset{j}\sum g(j)}
i∏j∏f(i)g(j)=i∏f(i)j∑g(j)
求和符号从指数跑下来就是连乘积符号,连乘积跑到指数上就是求和。
对这个式子做一下处理:
∏
i
=
1
n
∏
j
=
1
m
f
(
g
c
d
(
i
,
j
)
)
\overset{n}{\underset{i=1}\prod}\overset{m}{\underset{j=1}\prod}f\left(gcd(i,j)\right)
i=1∏nj=1∏mf(gcd(i,j))
首先发现无法可想,可暴力枚举一下
g
c
d
gcd
gcd:
=
∏
g
=
1
n
∏
i
=
1
n
∏
j
=
1
m
f
(
g
)
[
g
c
d
(
i
,
j
)
=
g
]
=\overset{n}{\underset{g=1}\prod}\overset{n}{\underset{i=1}\prod}\overset{m}{\underset{j=1}\prod}f\left(g\right)^{[gcd(i,j)=g]}
=g=1∏ni=1∏nj=1∏mf(g)[gcd(i,j)=g]
=
∏
g
=
1
n
∏
i
=
1
⌊
n
g
⌋
∏
j
=
1
⌊
m
g
⌋
f
(
g
)
[
g
c
d
(
i
,
j
)
=
1
]
=\overset{n}{\underset{g=1}\prod}\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\prod}\overset{\left\lfloor\frac m{g}\right\rfloor}{\underset{j=1}\prod}f\left(g\right)^{[gcd(i,j)=1]}
=g=1∏ni=1∏⌊gn⌋j=1∏⌊gm⌋f(g)[gcd(i,j)=1]
把
p
r
o
d
prod
prod式推到指数上:
=
∏
g
=
1
n
f
(
g
)
∑
i
=
1
⌊
n
g
⌋
∑
j
=
1
⌊
m
g
⌋
[
g
c
d
(
i
,
j
)
=
1
]
=\overset{n}{\underset{g=1}\prod}f\left(g\right)^{\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}\overset{\left\lfloor\frac m{g}\right\rfloor}{\underset{j=1}\sum}[gcd(i,j)=1]}
=g=1∏nf(g)i=1∑⌊gn⌋j=1∑⌊gm⌋[gcd(i,j)=1]
对指数做一波反演…:
=
∏
g
=
1
n
f
(
g
)
∑
d
=
1
n
μ
(
d
)
⌊
n
d
g
⌋
⌊
m
d
g
⌋
=\overset{n}{\underset{g=1}\prod}f\left(g\right)^{\overset{n}{\underset{d=1}\sum}\mu(d)\left\lfloor\frac n{dg}\right\rfloor\left\lfloor\frac m{dg}\right\rfloor}
=g=1∏nf(g)d=1∑nμ(d)⌊dgn⌋⌊dgm⌋
似乎无法可想了,所以再把和式拉下来:
=
∏
g
=
1
n
∏
d
=
1
n
f
(
g
)
μ
(
d
)
⌊
n
d
g
⌋
⌊
m
d
g
⌋
=\overset{n}{\underset{g=1}\prod}\overset{n}{\underset{d=1}\prod}f\left(g\right)^{\mu(d)\left\lfloor\frac n{dg}\right\rfloor\left\lfloor\frac m{dg}\right\rfloor}
=g=1∏nd=1∏nf(g)μ(d)⌊dgn⌋⌊dgm⌋
换元,令
T
=
d
p
T=dp
T=dp:
=
∏
T
=
1
n
∏
d
=
1
n
f
(
T
d
)
μ
(
d
)
[
d
∣
T
]
⌊
n
T
⌋
⌊
m
T
⌋
=\overset{n}{\underset{T=1}\prod}\overset{n}{\underset{d=1}\prod}f\left(\frac T d\right)^{\mu(d)[d|T]\left\lfloor\frac n{T}\right\rfloor\left\lfloor\frac m{T}\right\rfloor}
=T=1∏nd=1∏nf(dT)μ(d)[d∣T]⌊Tn⌋⌊Tm⌋
=
∏
T
=
1
n
(
∏
d
=
1
n
f
(
T
d
)
μ
(
d
)
[
d
∣
T
]
)
⌊
n
T
⌋
⌊
m
T
⌋
=\overset{n}{\underset{T=1}\prod}\left(\overset{n}{\underset{d=1}\prod}f\left(\frac T d\right)^{\mu(d)[d|T]}\right)^{\left\lfloor\frac n{T}\right\rfloor\left\lfloor\frac m{T}\right\rfloor}
=T=1∏n(d=1∏nf(dT)μ(d)[d∣T])⌊Tn⌋⌊Tm⌋
每个
d
d
d都只对它的倍数做贡献,括号内求一个区间积,括号外面整除分块,于是就做完了。
对
g
g
g换元当然也是可做的。
#include<iostream>
#include<vector>
using namespace std;
long long pow(long long a,long long b,long long c) {
long long ans=1;
while(b) {
if(b&1) ans=(ans*a)%c;
a=(a*a)%c;
b>>=1;
}
return ans;
}
const int N=1e6;
const long long M=1e9+7;
bool vis[N+5];
vector<int> sta;
int mu[N+5];
long long f[N+5],g[N+5],F[N+5];
void Euler() {
mu[1]=1;
for(int i=2; i<=N; i++) {
if(!vis[i])
sta.push_back(i),
mu[i]=-1;
for(auto&j:sta) {
if(i*j>N) break;
int m=i*j;
vis[m]=true;
if(i%j) mu[m]=-mu[i];
else break;
}
}
F[0]=F[1]=f[1]=g[1]=1;
for(int i=2; i<=N; i++)
f[i]=(f[i-1]+f[i-2])%M,
g[i]=pow(f[i],M-2,M),
F[i]=1;
for(int i=1; i<=N; i++)
for(int j=i; j<=N;j+=i)
if(mu[j/i])
F[j]=F[j]*(~mu[j/i]?f[i]:g[i])%M;
for(int i=1; i<=N; i++)
(F[i]*=F[i-1])%=M;
}
long long calc(int n,int m) {
long long ans=1;
long long t;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),
t=F[r]*pow(F[l-1],M-2,M)%M,
(ans*=pow(t,(long long)(n/l)*(m/l)%(M-1),M))%=M;//这里必须要括起来是因为整除
return ans;
}
int main() {
int n,m;
Euler();
// for(int i=1;i<=10;i++)
// cout<<F[i]<<' ';
int T;
cin>>T;
while(T--) {
cin>>n>>m;
if(n>m) swap(n,m);
cout<<calc(n,m)<<endl;
}
}
狄利克雷卷积,杜教筛
题目要求:
∑
i
=
1
n
∑
j
=
1
n
i
j
gcd
(
i
,
j
)
\overset{n}{\underset{i=1}\sum}\overset{n}{\underset{j=1}\sum}ij\gcd(i,j)
i=1∑nj=1∑nijgcd(i,j)
=
∑
g
=
1
n
g
∑
i
=
1
n
i
∑
j
=
1
n
j
[
g
c
d
(
i
,
j
)
=
g
]
=\overset{n}{\underset{g=1}\sum}g\overset{n}{\underset{i=1}\sum}i\overset{n}{\underset{j=1}\sum}j[gcd(i,j)=g]
=g=1∑ngi=1∑nij=1∑nj[gcd(i,j)=g]
=
∑
g
=
1
n
g
3
∑
i
=
1
⌊
n
g
⌋
i
∑
j
=
1
⌊
n
g
⌋
j
[
g
c
d
(
i
,
j
)
=
1
]
=\overset{n}{\underset{g=1}\sum}g^3\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}i\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{j=1}\sum}j[gcd(i,j)=1]
=g=1∑ng3i=1∑⌊gn⌋ij=1∑⌊gn⌋j[gcd(i,j)=1]
=
∑
d
=
1
n
μ
(
d
)
∑
g
=
1
n
g
3
∑
i
=
1
⌊
n
g
⌋
[
d
∣
i
]
i
∑
j
=
1
⌊
n
g
⌋
[
d
∣
j
]
j
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{n}{\underset{g=1}\sum}g^3\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{i=1}\sum}[d|i]i\overset{\left\lfloor\frac n{g}\right\rfloor}{\underset{j=1}\sum}[d|j]j
=d=1∑nμ(d)g=1∑ng3i=1∑⌊gn⌋[d∣i]ij=1∑⌊gn⌋[d∣j]j
=
∑
d
=
1
n
μ
(
d
)
∑
g
=
1
n
d
2
g
3
(
∑
i
=
1
⌊
n
d
g
⌋
i
)
2
=\overset{n}{\underset{d=1}\sum}\mu(d)\overset{n}{\underset{g=1}\sum}d^2g^3\left(\overset{\left\lfloor\frac n{dg}\right\rfloor}{\underset{i=1}\sum}i\right)^2
=d=1∑nμ(d)g=1∑nd2g3
i=1∑⌊dgn⌋i
2
换元,
T
=
d
g
T=dg
T=dg:
=
∑
T
=
1
n
∑
g
=
1
n
g
μ
(
T
g
)
[
g
∣
T
]
(
∑
i
=
1
⌊
n
T
⌋
i
)
2
=\overset{n}{\underset{T=1}\sum}\overset{n}{\underset{g=1}\sum}g\mu(\frac T g)[g|T]\left(\overset{\left\lfloor\frac n{T}\right\rfloor}{\underset{i=1}\sum}i\right)^2
=T=1∑ng=1∑ngμ(gT)[g∣T]
i=1∑⌊Tn⌋i
2
=
∑
T
=
1
n
∑
g
∣
T
μ
(
T
g
)
g
(
∑
i
=
1
⌊
n
T
⌋
i
)
2
=\overset{n}{\underset{T=1}\sum}\overset{}{\underset{g|T}\sum}\mu(\frac T g)g\left(\overset{\left\lfloor\frac n{T}\right\rfloor}{\underset{i=1}\sum}i\right)^2
=T=1∑ng∣T∑μ(gT)g
i=1∑⌊Tn⌋i
2
=
∑
T
=
1
n
(
∑
i
=
1
⌊
n
T
⌋
i
)
2
∑
g
∣
T
μ
(
T
g
)
g
=\overset{n}{\underset{T=1}\sum}\left(\overset{\left\lfloor\frac n{T}\right\rfloor}{\underset{i=1}\sum}i\right)^2\overset{}{\underset{g|T}\sum}\mu(\frac T g)g
=T=1∑n
i=1∑⌊Tn⌋i
2g∣T∑μ(gT)g
最后一个和式很明显是狄利克雷卷积:
=
∑
T
=
1
n
(
∑
i
=
1
⌊
n
T
⌋
i
)
2
φ
(
T
)
=\overset{n}{\underset{T=1}\sum}\left(\overset{\left\lfloor\frac n{T}\right\rfloor}{\underset{i=1}\sum}i\right)^2\varphi(T)
=T=1∑n
i=1∑⌊Tn⌋i
2φ(T)
杜教筛一波就做完了…这个题好像就必须用到小学奥数了,自己看着办吧。
#include<iostream>
#include<unordered_map>
#include<vector>
using namespace std;
const int N=1e7;
long long M;
long long phi[N+5];
bool vis[N+5];
vector<int> sta;
long long pow(long long a,long long b,long long c) {
long long ans=1;
while(b) {
if(b&1) ans=(ans*a)%c;
a=(a*a)%c;
b>>=1;
}
return ans;
}
void Euler() {
phi[1]=1;
for(int i=2;i<=N;i++) {
if(!vis[i])
phi[i]=i-1,
sta.push_back(i);
for(auto&j:sta) {
long long m=(long long)i*j;
if(m>N) break;
vis[m]=true;
if(i%j)
phi[m]=phi[i]*phi[j];
else {
phi[m]=phi[i]*j;
break;
}
}
}
for(long long i=1;i<=N;i++)
(phi[i]=phi[i]*i%M*i%M+phi[i-1])%=M;
}
long long inv;
long long f(long long n) {
n%=M;
return (n*(n+1)/2)%M*((n*(n+1)/2)%M)%M;
}
long long h(long long n) {
n%=M;
return n*(n+1)%M*(2*n+1)%M*inv%M;
}
//typedef long long LL;
//LL h(LL x){x%=M;return x*(x+1)%M*(2*x+1)%M*inv%M;}
//LL f(LL x){x%=M;return(x*(x+1)/2)%M*((x*(x+1)/2)%M)%M;}
unordered_map<long long,long long> H;
long long sum(long long n) {
if(n<=N) return phi[n];
if(H.count(n)) return H[n];
long long ans=f(n);
for(long long l=2,r;l<=n;l=r+1)
r=n/(n/l),
ans=(ans-(h(r)-h(l-1))%M*sum(n/l)%M+M)%M;
return H[n]=ans;
}
long long calc(long long n) {
long long ans=0;
for(long long l=1,r;l<=n;l=r+1)
r=n/(n/l),
ans=(ans+(sum(r)-sum(l-1))%M*f(n/l)%M+M)%M;
return ans;
}
signed main() {
long long n;
cin>>M>>n;
inv=pow(6,M-2,M);
// cout<<"***"<<inv<<endl;
Euler();
// cout<<M<<endl;
// cout<<n;
// cout<<sum(1e10)<<endl;
// for(int i=1;i<=10;i++)
// cout<<sum(i+1e6)<<' ';
// cout<<endl;
cout<<calc(n);
}
两种换元方式都可以,用欧拉反演这道题会简单许多,上一道题用不了欧拉反演的原因就是,它把 g c d gcd gcd搞到了分母上。
递推
其中
d
d
d函数就是整除函数
σ
0
σ_0
σ0.
首先注意到整除函数有如下性质:
那么我们可以替换一下题目中的式子:
这就是一个莫反的形式,一种做法是同时对三个艾弗森括号进行莫比乌斯反演,最后枚举三元环,可以通过此题。
另一种做法是只对一个艾弗森括号进行莫比乌斯反演,这是数论做法,也可以通过此题。
这里选择对
[
y
⊥
z
]
[y⊥z]
[y⊥z]做反演:
然后整理下式子:
然后在这里,一般来说是把
d
d
d拉到最左边,不过这一次是把
x
x
x拉到最左边,然后继续整理:
注意这里艾弗森括号可以拆开。
然后发现
[
x
⊥
d
]
[x⊥d]
[x⊥d]这一项与后面两个和式没关系了,可以拉出来:
注意到艾弗森括号只有
0
0
0或者
1
1
1两种取值,所以平方符号可以不写:
然后呢?似乎无法可想。
所以可以考虑设函数了:
这时候可以考虑数论分块了:
此时这两个函数如果暴算的话,我们枚举
x
x
x和
d
d
d,总复杂度就是O(n2logn),因为
f
f
f的
N
N
N取值是一个调和级数。
当然还有一个神仙的做法,其实这两个函数是可以递推的:
首先假设
k
k
k不含有任何相同质因子。也就是说
μ
(
k
)
≠
0
μ(k)≠0
μ(k)=0。
如果有
k
′
k'
k′含有相同质因子,那么显然有对应的
k
k
k只含有
k
′
k'
k′的每个质因子一次,有
f
(
N
,
k
)
=
f
(
N
,
k
′
)
,
g
(
N
,
k
)
=
g
(
N
,
k
′
)
f(N,k)=f(N,k'),g(N,k)=g(N,k')
f(N,k)=f(N,k′),g(N,k)=g(N,k′)。因为只需要判断是否含有相同质因子,就能判断是否互质,而不关乎相同质因子数目。
那么对于上面的
k
k
k,其实只需要做一个简单的容斥:
f
(
N
,
k
)
=
∑
i
=
1
N
⌊
N
i
⌋
[
i
⊥
k
]
f(N,k)=\overset{N}{\underset{i=1}\sum}\left\lfloor\frac Ni\right\rfloor[i\perp k]
f(N,k)=i=1∑N⌊iN⌋[i⊥k]
=
∑
i
=
1
N
⌊
N
i
⌋
[
i
⊥
k
p
]
−
∑
i
=
1
N
⌊
N
i
⌋
[
i
⊥
k
p
]
[
p
∣
i
]
=\overset{N}{\underset{i=1}\sum}\left\lfloor\frac Ni\right\rfloor[i\perp \frac kp]-\overset{N}{\underset{i=1}\sum}\left\lfloor\frac Ni\right\rfloor[i\perp \frac kp][p|i]
=i=1∑N⌊iN⌋[i⊥pk]−i=1∑N⌊iN⌋[i⊥pk][p∣i]
=
∑
i
=
1
N
⌊
N
i
⌋
[
i
⊥
k
p
]
−
∑
i
=
1
N
⌊
N
i
⌋
[
i
⊥
k
p
]
[
p
∣
i
]
=\overset{N}{\underset{i=1}\sum}\left\lfloor\frac Ni\right\rfloor[i\perp \frac kp]-\overset{N}{\underset{i=1}\sum}\left\lfloor\frac Ni\right\rfloor[i\perp \frac kp][p|i]
=i=1∑N⌊iN⌋[i⊥pk]−i=1∑N⌊iN⌋[i⊥pk][p∣i]
=
∑
i
=
1
N
⌊
N
i
⌋
[
i
⊥
k
p
]
−
∑
i
=
1
⌊
N
p
⌋
⌊
N
p
⋅
i
⌋
[
p
⋅
i
⊥
k
p
]
=\overset{N}{\underset{i=1}\sum}\left\lfloor\frac Ni\right\rfloor[i\perp \frac kp]-\overset{\left\lfloor\frac Np\right\rfloor}{\underset{i=1}\sum}\left\lfloor\frac N{p\cdot i}\right\rfloor[p\cdot i\perp \frac kp]
=i=1∑N⌊iN⌋[i⊥pk]−i=1∑⌊pN⌋⌊p⋅iN⌋[p⋅i⊥pk]
注意到一定有
p
⊥
k
p
p\perp\frac kp
p⊥pk,因此可以去掉:
=
∑
i
=
1
N
⌊
N
i
⌋
[
i
⊥
k
p
]
−
∑
i
=
1
⌊
N
p
⌋
⌊
N
p
⋅
i
⌋
[
i
⊥
k
p
]
=\overset{N}{\underset{i=1}\sum}\left\lfloor\frac Ni\right\rfloor[i\perp \frac kp]-\overset{\left\lfloor\frac Np\right\rfloor}{\underset{i=1}\sum}\left\lfloor\frac N{p\cdot i}\right\rfloor[i\perp \frac kp]
=i=1∑N⌊iN⌋[i⊥pk]−i=1∑⌊pN⌋⌊p⋅iN⌋[i⊥pk]
=
f
(
N
,
k
p
)
−
f
(
⌊
N
p
⌋
,
k
p
)
=f\left(N,\frac kp\right)-f\left(\left\lfloor\frac Np\right\rfloor,\frac kp\right)
=f(N,pk)−f(⌊pN⌋,pk)
对于
g
g
g,也有类似的容斥:
f
(
N
,
k
)
=
∑
i
=
1
N
μ
(
i
)
[
i
⊥
k
]
f(N,k)=\overset{N}{\underset{i=1}\sum}\mu(i)[i\perp k]
f(N,k)=i=1∑Nμ(i)[i⊥k]
=
∑
i
=
1
N
μ
(
i
)
[
i
⊥
k
p
]
−
∑
i
=
1
N
μ
(
i
)
[
i
⊥
k
p
]
[
p
∣
i
]
=\overset{N}{\underset{i=1}\sum}\mu(i)[i\perp \frac kp]-\overset{N}{\underset{i=1}\sum}\mu(i)[i\perp \frac kp][p|i]
=i=1∑Nμ(i)[i⊥pk]−i=1∑Nμ(i)[i⊥pk][p∣i]
=
∑
i
=
1
N
μ
(
i
)
[
i
⊥
k
p
]
−
∑
i
=
1
⌊
N
p
⌋
μ
(
p
⋅
i
)
[
i
⊥
k
p
]
=\overset{N}{\underset{i=1}\sum}\mu(i)[i\perp \frac kp]-\overset{\left\lfloor\frac Np\right\rfloor}{\underset{i=1}\sum}\mu(p\cdot i)[i\perp \frac kp]
=i=1∑Nμ(i)[i⊥pk]−i=1∑⌊pN⌋μ(p⋅i)[i⊥pk]
若
p
⊥
i
p\perp i
p⊥i,则
μ
(
p
⋅
i
)
=
−
μ
(
i
)
\mu(p\cdot i)=-\mu(i)
μ(p⋅i)=−μ(i),否则
μ
(
p
⋅
i
)
=
0
\mu(p\cdot i)=0
μ(p⋅i)=0:
=
∑
i
=
1
N
μ
(
i
)
[
i
⊥
k
p
]
+
∑
i
=
1
⌊
N
p
⌋
μ
(
i
)
[
i
⊥
p
]
[
i
⊥
k
p
]
=\overset{N}{\underset{i=1}\sum}\mu(i)[i\perp \frac kp]+\overset{\left\lfloor\frac Np\right\rfloor}{\underset{i=1}\sum}\mu(i)[i\perp p][i\perp \frac kp]
=i=1∑Nμ(i)[i⊥pk]+i=1∑⌊pN⌋μ(i)[i⊥p][i⊥pk]
i
⊥
p
i\perp p
i⊥p与
i
⊥
k
p
i\perp \frac kp
i⊥pk显然可以合起来:
=
∑
i
=
1
N
μ
(
i
)
[
i
⊥
k
p
]
+
∑
i
=
1
⌊
N
p
⌋
μ
(
i
)
[
i
⊥
k
]
=\overset{N}{\underset{i=1}\sum}\mu(i)[i\perp \frac kp]+\overset{\left\lfloor\frac Np\right\rfloor}{\underset{i=1}\sum}\mu(i)[i\perp k]
=i=1∑Nμ(i)[i⊥pk]+i=1∑⌊pN⌋μ(i)[i⊥k]
=
g
(
N
,
k
p
)
+
g
(
⌊
N
p
⌋
,
k
)
=g(N,\frac kp)+g\left(\left\lfloor\frac Np\right\rfloor,k\right)
=g(N,pk)+g(⌊pN⌋,k)
我们再回顾一下递推式:
然后是初值:
所有的
f
(
x
,
1
)
f(x,1)
f(x,1)的值,可以数论分块,O(
N
N
N^1.5)求出。
当然还有一种做法:
这样可以O(
N
N
N)求出所有的
f
(
x
,
1
)
f(x,1)
f(x,1)。
这个式子可以证明一下,从左到右证明较难,而从右到左证明较易:
因此:
那么我们分析一下复杂度:
f
,
g
f,g
f,g的第一维,取遍
N
N
N的取整点,复杂度为
O
(
N
)
O(\sqrt N)
O(N),第二维为
μ
(
i
)
≠
0
\mu(i)\neq 0
μ(i)=0的数量,
N
=
100000
N=100000
N=100000时爆搜一波有
60000
60000
60000个,
O
(
N
)
O(N)
O(N)。
然后就是调用 g ( l − 1 ) , g ( r ) g(l-1),g(r) g(l−1),g(r)的问题,由于有 O ( n ) O(\sqrt n) O(n)个 l − 1 , r l-1,r l−1,r,每个 l − 1 , r l-1,r l−1,r都会对应 O ( r ) O(\sqrt r) O(r)个取整点,会不会使得复杂度增加为 O ( n ) O(n) O(n)的问题。答案是不会,我们只需要讨论 r r r,因为 l − 1 l-1 l−1肯定就是上一个块的 r r r。我们注意到 r = n / ( n / l ) r=n/(n/l) r=n/(n/l), r r r本身就是 n n n的取整点。
因此写了一波二维unordered_map优化整除分块,被卡成30。考虑递推,先对质因子个数少的 x x x统计答案,再根据质因子个数少的 f , g f,g f,g函数值计算出质因子个数多的 x x x的 f , g f,g f,g函数值,然后再计算答案:
#include<iostream>
#include<vector>
#include<unordered_map>
using namespace std;
const int N=1e5;
const long long M=1e9+7;
vector<int> sta;
bool vis[N+5];
int mu[N+5],LPF[N+5],con[N+5],s0[N+5],LPFx[N+5],cnt[N+5];
//莫比乌斯函数,最小质因子,本质相同素数,约数个数,含有的最小质因子的次数,本质不同质因子个数+1
int now;
long long f[N+5][20],g[N+5][20];
void Euler() {
mu[1]=1;s0[1]=1;LPF[1]=1;LPFx[1]=1;con[1]=1;cnt[1]=1;
for(int i=2; i<=N; i++) {
if(!vis[i])
mu[i]=-1,LPF[i]=i,con[i]=i,LPFx[i]=1,s0[i]=2,cnt[i]=2,sta.push_back(i);
for(auto&j:sta) {
int x=i*j;
if(x>N) break;
vis[x]=true;
LPF[x]=j;
if(i%j)
mu[x]=-mu[i],con[x]=con[i]*con[j],s0[x]=s0[i]*s0[j],LPFx[x]=1,cnt[x]=cnt[i]+1;
else {
con[x]=con[i],LPFx[x]=LPFx[i]+1,cnt[x]=cnt[i],s0[x]=s0[i]/LPFx[x]*(LPFx[x]+1);
break;
}
}
}
for(int i=1; i<=N; i++) f[i][0]=g[i][0]=1,f[i][1]=f[i-1][1]+s0[i],g[i][1]=g[i-1][1]+mu[i];
}
inline __attribute__((always_inline)) __uint128_t MOD(__uint128_t&x) {
return x>=1e19?x%=M:M;
}
int n,m,K;
__uint128_t d[N+5];
long long calc(int x,int t) {
int k=cnt[x];
__uint128_t sum=0;
for(int l=1,r;l<=m;l=r+1)
r=min(m/(m/l),K/(K/l)),
MOD(sum+=(g[r][k]-g[l-1][k])*f[m/l][k]*f[K/l][k]);
MOD(sum*=d[x]);
for(int i=t;i<(int)sta.size()&&(long long)x*sta[i]<=n;i++) {
for(int l=1,r;l<=m;l=r+1)
r=min(m/(m/l),K/(K/l)),
g[r][k+1]=g[r][k]+g[r/sta[i]][k+1],//g只计算右端点就可以了,左端点-1就是上一块的右端点
f[m/l][k+1]=f[m/l][k]-f[m/l/sta[i]][k],
f[K/l][k+1]=f[K/l][k]-f[K/l/sta[i]][k];
MOD(sum+=calc(x*sta[i],i+1));
}
return sum%M;
}
int main() {
Euler();
int T;
cin>>T;
while(T--) {
cin>>n>>m>>K;
if(n>K) swap(n,K);
if(n>m) swap(n,m);
if(m>K) swap(m,K);
for(int i=1;i<=n;i++) d[i]=0;
for(int i=1;i<=n;i++) d[con[i]]+=n/i;
cout<<calc(1,0)<<endl;
}
return 0;
}
后记
于是皆大欢喜。