「笔记」莫比乌斯反演
Updated on 2020.8.6
巨幅更新,对积性函数和狄利克雷卷积部分进行重构。
新增对一类特殊求和式的讲解。
Updated on 2020.9.9
添加了几个例题。
前置知识
小碎骨
艾佛森括号 \([P] = \begin{cases} 1 &\text{If P is true}\\ 0 &\text{Otherwise} \end{cases}\)
此处 \(P\) 是一个可真可假的命题。
引理1
证明
数论分块
内容独立了出来,详细内容见 数论分块 - Luckyblock
对于一类含有\(\left\lfloor\frac{n}{i}\right\rfloor\)的求和式 (\(n\) 为常数),由于\(\left\lfloor\frac{n}{i}\right\rfloor\)单调不增,故存在多个区间\([l,r]\), 使得\(\left\lfloor\frac{n}{i}\right\rfloor = \left\lfloor\frac{n}{j}\right\rfloor(i,j\in [l,r])\) 成立。
对于任意一个\(i\),最大的满足上式的 \(j=\left\lfloor{\dfrac{n}{\left\lfloor{\dfrac{n}{i}}\right\rfloor}}\right\rfloor\)
积性函数
定义
若 \(\gcd(x,y) = 1\) 且 \(f(xy)=f(x)f(y)\),则 \(f(n)\) 为积性函数。
性质
若 \(f(x)\),\(g(x)\)均为积性函数,则以下函数也为积性函数:
常见积性函数
- 单位函数 \(e(n) = [n = 1]\)。
- 幂函数 \(\operatorname{Id}_{k}(n) = n^k\), \(\operatorname{id}_1(n)\) 通常简记为\(\operatorname{id}(n)\)。
- 常数函数 \(1(n) = 1\)。
- 因数个数 \(\operatorname{d}(n) = \sum\limits_{d\mid n} 1\)。
- 除数函数 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)。
\(k=1\) 时为因数和函数,通常简记为 \(\sigma(n)\)。
\(k=0\) 时为因数个数函数 \(\sigma_{0}(n)\)。 - 欧拉函数 \(\varphi(n) = \sum\limits_{i=1}^{n} [\gcd(i,n) = 1]\)。
- 莫比乌斯函数 \(\mu(n) = \begin{cases}1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{为}\ n\ \text{的本质不同质因子个数} \end{cases}\)
不是很懂上面写的什么玩意?
不用深究,有个印象继续往下看就好。
莫比乌斯函数
定义
\(\mu\) 为莫比乌斯函数,定义为
解释
令 \(n = \prod\limits_{i=1}^{k} p_{i}^{c_i}\),其中\(p_i\)为质因子,\(c_i\ge 1\)。
- \(n=1\)时,\(\mu (n) = 1\)。
- \(n\not ={1}\)时 ,
- \(\exist i\in [1,k], c_i > 1\) 时,\(\mu (n) = 0\)。
当某质因子出现次数大于\(1\)时,\(\mu (n) = 0\)。 - \(\forall i\in [1,k], c_i = 1\) 时,\(\mu (n) = (-1)^k\)。
当每个质因子只出现一次时,即\(n = \prod\limits_{i=1}^{k}p_i\),\(\{p_i\}\)中元素唯一。
\(\mu (n) = (-1)^k\),此处\(k\)为质因子的种类数。
- \(\exist i\in [1,k], c_i > 1\) 时,\(\mu (n) = 0\)。
性质
莫比乌斯函数是积性函数,且具有以下性质
证明,设 \(n = \prod\limits_{i=1}^{k}{p_i^{c_i}}, n' = \prod\limits_{i=1}^{k}{p_i}\)。
- 根据莫比乌斯函数定义,则有:\(\sum\limits_{d\mid n}{\mu(d)} = \sum\limits_{d\mid n'}{\mu(d)}\)。
- 若 \(n'\) 的某因子 \(d\),有 \(\mu (d) = (-1)^i\),则它由 \(i\) 个 本质不同的质因子组成。
由于质因子总数为 \(k\),满足上式的因子数为 \(C_{k}^{i}\)。 - 对于原求和式,转为枚举 \(\mu(d)\) 的值。
则 \(\sum\limits_{d\mid n'}{\mu(d)} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i\times 1^{k-i}}\)
根据二项式定理,上式 \(= (1+(-1))^k\),
易知该式在 \(k=0\),即 \(n=1\) 时为 \(1\),否则为 \(0\)。
反演常用结论
一个反演常用结论:
证明 1:
设 \(n = \gcd(i,j)\),则右\(= \sum\limits_{d\mid n} {\mu (d)} = [n = 1] = [\gcd(i,j) = 1]=\) 左。
证明 2:
暴力展开:\([\gcd(i,j) = 1] = e(\gcd(i,j)) = \sum\limits_{d\mid \gcd(i,j)}\mu(d)\)
线性筛求莫比乌斯函数
\(\mu\) 为积性函数,因此可以线性筛莫比乌斯函数。
int pnum, mu[kMaxn], p[kMaxn];
bool vis[kMaxn];
void Euler(int lim_) {
vis[1] = true, mu[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = - mu[i];
}
}
}
狄利克雷(Dirichlet)卷积
建议阅读 算法学习笔记(35): 狄利克雷卷积 By: Pecco。
定义两个数论函数 \(f,g\) 的狄利克雷卷积为
性质
- 显然满足 交换律,结合律,分配律。
- \(f \ast g = g \ast f\)
- \((f \ast g) \ast h = f\ast (g\ast h)\)
- \(f\ast (g+h) = f\ast g + f\ast h\)
- \(e\) 为狄利克雷卷积的单位元,有\((f\ast e)(n) = f(n)\)。
- 若 \(f,g\) 为积性函数,则 \(f\ast g\) 为积性函数。
关于单位元 \(e\)
有:
证明
- 对于\([\dfrac{n}{d} = 1]\),当且仅当 \(\dfrac{n}{d}=1\),即 \(d=n\) 时为 \(1\),否则为\(0\)。
- 则当 \(d=n\) 时,\(f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\)。
当 \(d\not ={n}\) 时,\(f(d)\left[\dfrac{n}{d} = 1\right] = 0\)。
综上,\((f\ast e)(n) = \sum\limits_{d\mid n} f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\),满足单位元的性质。
则 \(e = \mu \ast 1\) 成立。
除数函数与幂函数
幂函数 \(\operatorname{Id}_{k}(n) = n^k\)。
除数函数 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)。
显然有:
当 \(k=0\) 时,\(\operatorname{Id_0}=1\),\(\sigma_0\) 为因数个数函数,有:
欧拉函数与恒等函数
对于一式,当 \(n=p^m\) 时(\(p\) 为质数),有:
\(p^i\) 的因子有 \(p^{i-1}\) 个,为 \(1\sim p^{i-1}\),故 \(\varphi(p^i) = p^i-p^{i-1}\)。
又 \((\varphi \ast 1)(n)\) 为积性函数,则对于任意正整数 \(n\),有:
得证。
对于 2 式,在 1 式基础上两侧同时 \(\ast \mu\) 即得。
左侧变为 \(\varphi \ast 1\ast \mu = \varphi \ast e = \varphi\)。
计算
考虑枚举 \(n\) 的因子,将贡献累加即可。
显然可以使用埃氏筛筛出所有前缀狄利克雷卷积,复杂度 \(O(nk\log n)\),其中 \(k\) 是计算函数值的复杂度。
莫比乌斯反演
反演是个啥?反演。
公式
设\(f(n),g(n)\)为两个数论函数。
如果有
那么有
证明
方法一:运用卷积。
原问题为:已知 \(f = g\ast 1\),证明 \(g = f\ast \mu\)。
易知如下转化:\(f\ast \mu = g\ast 1 \ast \mu \Longrightarrow f\ast \mu = g\ast e = g\)。
方法二:对原式进行数论变换。
-
用 \(\sum\limits_{d\mid n}g(d)\) 替换\(f\left(\dfrac{n}{d}\right)\)。
\[\sum_{d\mid n}{\mu(d)\sum_{k\mid \frac{n}{d}}g(k)} \] -
变换求和顺序。
\[\sum_{k\mid n}g(k)\sum_{d\mid \frac{n}{k}}{\mu(d)} \] -
依据 \(\sum\limits_{d\mid n}{\mu(d)} = [n=1]\),仅当 \(\dfrac{n}{k} = 1\) 时,\(\sum\limits_{d\mid \frac{n}{k}}{\mu(d)} = 1\),否则为 \(0\)。
此时\(k=n\),故上式等价于 \(\sum\limits_{k\mid n} {[n=k]\cdot g(k)} = g(n)\)。
举例
从 狄利克雷(Dirichlet)卷积 部分可以知道一些积性函数的关系。
尝试对它们进行反演:
关于一类求和式
一般套路
考虑构造出函数 \(g\),满足下列形式:
则求和式变为:
考虑算术基本定理,发现若满足 \(d\mid \gcd (i,j)\),则 \(d\mid i\) 且 \(d\mid j\) 成立。
考虑 \(g(d)\) 在何时做出贡献,调整枚举顺序:
\(\sum\limits_{i=1}^{n}[d\mid i]\) 等价于 \(1\sim n\) 中 \(d\) 的倍数的个数,则上式等价于:
数论分块求解即可。
例 1
发现此题的 \(f\) 等价于 \(\operatorname{Id}\),则上式等价于:
例 2
感悟
卷点什么东西,把 \(g\) 卷出来。
\(g\) 不一定是特殊意义的函数。
例题
[HAOI2011] Problem b
\(n\) 组询问,每次给定参数 \(a,b,c,d,k\),求:
\[\sum\limits_{x=a}^{b}\sum\limits_{y=c}^{d}[\gcd(x,y) = k] \]\(1\le n,k,a,b,c,d\le 5\times 10^4\),\(a\le b,c\le d\)。
令 \(f(n,m) = \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = k]\)。
根据容斥原理,则原式等价于 \(f(b,d) - f(a-1,d) - f(b,d-1) + f(a-1,d-1)\)。
\(f\) 变成了上述一类求和式的形式,考虑化简 \(f\)。
易知原式等价于
代入反演常用结论 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),上式化为
变换求和顺序,先枚举\(d\mid \gcd(i,j)\),可得
对于上式后半的解释:当\(d\mid i\) 且 \(d\mid j\) 时,\(d\mid \gcd(i,j)\)。
易知\(1\sim \left\lfloor\dfrac{n}{k}\right\rfloor\) 中 \(d\) 的倍数有 \(\left\lfloor\dfrac{\left\lfloor\dfrac{n}{k}\right\rfloor}{d}\right\rfloor = \left\lfloor\dfrac{n}{kd}\right\rfloor\) 个(由引理 1 可知),原式变为
预处理 \(\mu\) 后,显然可以数论分块求解,复杂度\(O(n + T\sqrt{n})\)。
代码
//知识点:莫比乌斯反演
/*
//By:Luckyblock
*/
#include <cstdio>
#include <cctype>
#include <algorithm>
#define ll long long
const int MARX = 6e4 + 10;
//=============================================================
int N, a, b, c, d, k;
int cnt, Mobius[MARX], Prime[MARX], Sum[MARX];
bool vis[MARX];
//=============================================================
inline int read()
{
int f = 1, w = 0; char ch = getchar();
for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = -1;
for(; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void GetMobius()
{
Mobius[1] = 1;
int MAX = MARX - 10;
for(int i = 2; i <= MAX; i ++)
{
if(! vis[i]) Mobius[i] = - 1, Prime[++ cnt] = i;
for(int j = 1; j <= cnt && i * Prime[j] <= MAX; j ++)
{
vis[i * Prime[j]] = true;
if(i % Prime[j] == 0) break;
Mobius[i * Prime[j]] = - Mobius[i];
}
}
for(int i = 1; i <= MAX; i ++) Sum[i] = Sum[i - 1] + Mobius[i];
}
ll Calc(int x, int y)
{
ll ans = 0ll; int max = std :: min(x, y);
for(int l = 1, r; l <= max; l = r + 1)
r = std :: min(x / (x / l), y / (y / l)),
ans += (1ll * x / (1ll * l * k)) * (1ll * y / (1ll * l * k)) * (Sum[r] - Sum[l - 1]);
return ans;
}
ll Solve()
{
a = read(), b = read(), c = read(), d = read(), k = read();
return Calc(b, d) - Calc(b, c - 1) - Calc(a - 1, d) + Calc(a - 1, c - 1);
}
//=============================================================
int main()
{
N = read(); GetMobius();
while(N --) printf("%lld\n", Solve());
return 0;
}
[国家集训队]Crash的数字表格
给定 \(n,m\) , 求:
\[\sum_{i=1}^n\sum_{j=1}^{m} \operatorname{lcm}(i,j)\bmod 20101009 \]\(1\le n,m\le 10^7\)。
易知原式等价于:
考虑枚举 \(\gcd(i,j)\),设枚举量为 \(d\)。
则 \(d=\gcd(i,j)\) 的充要条件是满足 \(d|i, d|j\) 且 \(\gcd(\dfrac{i}{d},\dfrac{j}{d}) = 1\),则原式等价于:
先枚举 \(d\),则原式等价于:
这个 \(d\) 很烦人,把 \(i,j\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\),\(\frac{j}{d}\)。
消去 \(d\mid i\),\(d\mid j\) 的限定条件,则原式等价于:
单独考虑后半部分,设 \(f(x,y) = \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\)。
发现 \(f(x,y)\) 的形式与 [HAOI2011] Problem b 中的式子类似,代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\) 套路一波:
前半部分 \(\sum\limits_{d=1}\mu(d) d^2\),可以考虑筛出 \(\mu(d)\) 后求前缀和。
后半部分是等差数列乘等差数列的形式,设 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\),\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算:
则对于 \(f(x,y)\),有:
数论分块求解即可。
再看回原式,原式等价于:
又是一个可以数论分块求解的形式。
线性筛预处理后 数论分块套数论分块,复杂度 \(O(n + m)\),瓶颈是线性筛。
一些注意的点
处理时会出现求平方的运算,需要特别注意取模问题,ll 都会爆,被坑惨了。
在预处理前缀和的这个地方:
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
注意先对平方取模,在乘上 \(\mu\),否则就会爆掉。
以及可以仅令 \(mu + mod\)。
以及这个地方:
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
平方计算,注意随时取模。
代码
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const ll kMod = 20101009;
const int kMaxn = 1e7 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], sum[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true, mu[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = - mu[i];
}
}
sum[1] = 1ll;
for (int i = 1; i <= lim_; ++ i) {
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
}
}
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
int f(int n_, int m_) {
int lim = std :: min(n_, m_), ret = 0;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
}
return ret;
}
int Sum(ll l, ll r) {
return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
}
//=============================================================
int main() {
int n = read(), m = read();
int lim = std :: min(n, m), ans = 0;
Euler(lim);
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
}
printf("%d", ans);
return 0;
}
/*
7718820 8445880
*/
SP5971 LCMSUM - LCM Sum
\(T\) 次询问,每次询问给定 \(n\),求:
\[\sum_{i=1}^{n}\operatorname{lcm}(i,n) \]\(1<T\le 3\times 10^5\),\(1\le n\le 10^6\)。
法一:无脑暴力
先拆 \(\operatorname{lcm}\),原式等价于:
套路的枚举 \(\gcd(i,n)\),调换枚举顺序,原式等价于:
把 \(i,n\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\),消去整除的条件,原式等价于:
代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),原式等价于:
值得注意的是 \(t\) 的上界为 \(\frac{n}{d}\),\(dt\le n\)。
调换枚举顺序,先枚举 \(t\),原式等价于:
套路地消去整除的条件,把 \(i\) 中的 \(t\) 提出来,原式等价于:
对于最后的一个求和项,设 \(g(x) = \sum\limits_{i=1}^{x}i = \frac{x(x+1)}{2}\),显然可以 \(O(1)\) 求解,原式等价于:
考虑枚举 \(T = dt\),显然 \(T\le n\)。
\(\mu(t)t\) 与 \(d\) 无关,可以直接考虑枚举 \(t|T\),原式等价于:
前半块是一个数论分块的形式,可以 \(O(\sqrt{n})\) 求解。
考虑后半块,设 \(f(n)=\sum\limits_{d|n}\mu(d)d\),发现它是一个积性函数,可以线性筛筛出,具体地:
其中 \(p\) 为 \(n\) 的最小质因子。
此时已经可以线性筛 + 数论分块求解,复杂度 \(O(n+T\sqrt{n})\),比较菜鸡,时限 500ms 过不了。
考虑筛出 \(f\) 后再用埃氏筛预处理 \(\sum\limits_{T=1}^{n}g(\left\lfloor\dfrac{n}{T}\right\rfloor)f(T)\),输出时乘上 \(n\),复杂度变为 \(O(n\log^2 n + n)\)。
法二:
同样先拆 \(\operatorname{lcm}\),枚举 \(\gcd(i,n)\),调换枚举顺序,原式等价于:
把 \(i,n\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\),消去整除的条件,原式等价于:
调整枚举对象,上式等价于:
考虑 \(\sum\limits_{i=1}^{d}[\gcd(i,d) = 1]i\) 的实际意义,表示 \([1,d]\) 中与 \(d\) 互质的数的和。
当 \(d>1\) 时,与 \(d\) 互质的数总是成对存在,即若 \(\gcd(i,d)=1\) 成立,则 \(\gcd(d-i,d)=1\) 成立。
每对这样的数的和为 \(d\),共有 \(\frac{\varphi(d)}{2}\) 对这样的数。
则原式等价于:
可以直接预处理答案。
预处理时先线性筛出 \(\varphi\),再埃氏筛枚举 \(i\) 的倍数,令它们的答案加上 \(\frac{\varphi(i)i}{2}\),最后输出时乘上 \(n\)。
复杂度 \(O(n\log^2 n + T)\)。
法二代码
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 1e6;
//=============================================================
ll phi[kMaxn + 10], ans[kMaxn + 10];
int pnum, p[kMaxn + 10];
bool flag[kMaxn + 10];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void GetMax(int &fir, int sec) {
if (sec > fir) fir = sec;
}
void GetMin(int &fir, int sec) {
if (sec < fir) fir = sec;
}
void GetPrime() {
phi[1] = 1, flag[1] = true; //注意初始化
for (int i = 2; i <= kMaxn; ++ i) {
if (! flag[i]) {
p[++ pnum] = i;
phi[i] = i - 1ll;
}
for (int j = 1; j <= pnum && i * p[j] <= kMaxn; ++ j) {
flag[i * p[j]] = true;
if (i % p[j]) {
phi[i * p[j]] = phi[i] * phi[p[j]];
} else {
phi[i * p[j]] = phi[i] * p[j];
break;
}
}
}
for (int i = 1; i <= kMaxn; ++ i) {
for (int j = 1; i * j <= kMaxn; ++ j) {
ans[i * j] += (i == 1 ? 1 : 1ll * phi[i] * i / 2);
}
}
}
//=============================================================
int main() {
GetPrime();
int T = read();
while (T --) {
int n = read();
printf("%lld\n", 1ll * ans[n] * n);
}
return 0;
}
[SDOI2015]约数个数和
\(T\) 次询问,每次询问给定 \(n,m\)。
定义 >\(\operatorname{d}(i)\) 为 \(i\) 的约数个数,求:\[\sum_{i=1}^{n}\sum_{j=1}^m\operatorname{d}(ij) \]\(1<T,n\le 5\times 10^4\)。
一个结论:
证明
先考虑 \(i = p^a\),\(j=p^b(p\in \mathrm{primes})\) 的情况,有:
对于等式左侧,\(p^{a+b}\) 的约数个数为 \(a+b+1\)。
对于等式右侧,在保证 \(\gcd(x,y)=1\) 成立的情况下,有贡献的数对 \((x,y)\) 只能是下列三种形式:
- \(x>0,y-0\),\(x\) 有 \(a\) 种取值方法。
- \(x=0,y>0\),\(y\) 有 \(b\) 种取值方法。
- \(x=0,y=0\)。
则等式右侧贡献次数为 \(a+b+1\) 次,等于 \(p^{a+b}\) 的约数个数。
则当 \(i = p^a\),\(j=p^b(p\in \mathrm{primes})\) 时等式成立。
又不同质因子间相互独立,上述情况可拓展到一般的情况。
对 \(\operatorname{d}(i,j)\) 进行化简,代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),原式等价于:
调换枚举顺序,先枚举 \(d\),原式等价于:
把各项中的 \(d\) 提出来,消去整除的条件,原式等价于:
将 \(\operatorname{d}(ij)\) 代回原式,原式等价于:
调换枚举顺序,先枚举 \(d\),原式等价于:
把 \(i,j\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}, \frac{j}{d}\),消去整除的条件,原式等价于:
考虑预处理 \(S(x) = \sum\limits_{i=1}^{x}\operatorname{d}(i)\),则原式等价于:
线性筛预处理 \(\mu,\operatorname{d}\),数论分块求解即可,复杂度 \(O(n+T\sqrt{n})\)。
代码
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 5e4 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 为最小质因子的次数
ll summu[kMaxn], sumd[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true;
mu[1] = d[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
num[i] = 1;
d[i] = 2;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
num[i * p[j]] = num[i] + 1;
d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
break;
}
mu[i * p[j]] = - mu[i];
num[i * p[j]] = 1;
d[i * p[j]] = 2ll * d[i]; //
}
}
for (int i = 1; i <= lim_; ++ i) {
summu[i] = summu[i - 1] + mu[i];
sumd[i] = sumd[i - 1] + d[i];
}
}
//=============================================================
int main() {
Euler(kMaxn - 10);
int T = read();
while (T --) {
int n = read(), m = read(), lim = std :: min(n, m);
ll ans = 0ll;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
}
printf("%lld\n", ans);
}
return 0;
}
/*
in
1
32 43
*/
/*
out
15420
*/
P3768 简单的数学题
给定参数 \(n\)、\(p\),求:
\[\left(\sum_{i=1}^n\sum_{j=1}^ni\cdot j\cdot \gcd(i,j)\right)\bmod p \]\(n\leq10^{10}\),\(5\times10^8\leq p\leq1.1\times10^9\),\(p\in \mathrm{primes}\)。
时限 4s。
无脑套路暴力。
考虑先枚举 \(\gcd(i,j)\),原式等价于:
提出 \(d\),变为枚举 \(\frac{i}{d}\) 与 \(\frac{j}{d}\),消去整除的条件,原式等价于:
代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),原式等价于:
值得注意的是 \(t\) 的上界为 \(\frac{n}{d}\),\(dt\le n\)。
调换枚举顺序,先枚举 \(t\),原式等价于:
和上面一样,提出 \(t\),套路地消去整除的条件,原式等价于:
发现后面两个求和是等差数列乘等差数列的形式。
设 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\),\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算:
代入原式,原式等价于:
考虑枚举 \(T = dt\),显然 \(T\le n\)。
再考虑枚举 \(d|T\),即可得到 \(t = \frac{T}{d}\),原式等价于:
对于后面这一坨,用 \(\sum\limits_{d|T}d\cdot \mu{\left(\frac{T}{d}\right)} = \operatorname{Id} \ast \mu(T)= \varphi(T)\) 反演,则原式等价于:
后半块可以数论分块,考虑前半块。
发现前半段即为 \(\operatorname{Id}^2(T)\varphi(T)\),又是前缀和形式,考虑杜教筛。
有:
考虑找到一个函数 \(g\),构造函数 \(h = f\ast g\) 使其便于求值,有:
看到同时存在 \(d\) 和 \(\frac{n}{d}\),考虑把 \(d^2\) 消去。
令 \(g = \operatorname{Id}^2\),有:
又 \(\varphi \ast 1 = \operatorname{Id}\),则有:
找到了合适的 \(g\),套杜教筛的公式。
前一项是自然数的立方和,有 \(\sum\limits_{i=1}^{n} i^3 = (\frac{n(n+1)}{2})^2\)。证明详见:自然数前n项平方和、立方和公式及证明 - 百度文库。
后一项直接等差数列求和 + 数论分块求解即可。
「SDOI2017」数字表格
设 \(f_{i}\) 表示斐波那契数列的第 \(i\) 项。
\(T\) 组数据,每次给定参数 \(n,m\),求:\[\prod_{i=1}^{n}\prod_{j=1}^{m}f_{\gcd(i,j)} \pmod {10^9 + 7} \]\(1\le T\le 10^3\),\(1\le n,m\le 10^6\)。
5S,256MB。
以下钦定 \(n\ge m\)。
大力化式子,先套路地枚举 \(\gcd(i,j)\),用初中知识把两个 \(\prod\) 化到指数位置,原式等于:
分母套路一波,有:
代回原式,原式等于:
考虑再暴力拆一波,原式等于:
做不动了,但发现变量仅有 \(k,d,kd\),考虑更换枚举对象改为枚举 \(t = kd\) 与 \(d\),则原式等于:
枚举对象变成了约数形式。从后面的式子推前面的式子是比较显然的,可以认为这种枚举 \(t=kd\) 的形式是一种逆向操作。
设:
\(g(t)\) 可以用类似埃氏筛的方法 \(O(n\log ^2 n)\) 地预处理出来。再把 \(g\) 代回原式,原式等于:
可以考虑预处理 \(g(t)\) 的前缀积,数论分块枚举指数求解即可。
总时间复杂度 \(O(n\log ^2 n + T\sqrt n)\),轻微卡常可以过。
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define LL long long
const LL mod = 1e9 + 7;
const int kN = 1e6;
//=============================================================
LL n, m, ans;
int p_num, p[kN + 10];
bool vis[kN + 10];
LL mu[kN + 10], f[kN + 10], g[kN + 10], prod[kN + 10];
LL invf[kN + 10], invp[kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) {
w = (w << 3) + (w << 1) + (ch ^ '0');
}
return f * w;
}
void Chkmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
x_ %= mod;
LL ret = 1;
for (; y_; y_ >>= 1ll) {
if (y_ & 1) ret = ret * x_ % mod;
x_ = x_ * x_ % mod;
}
return ret;
}
void Euler() {
vis[1] = true, mu[1] = 1;
for (int i = 2; i <= kN; ++ i) {
if (! vis[i]) {
p[++ p_num] = i;
mu[i] = -1;
}
for (int j = 1; j <= p_num && i * p[j] <= kN; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
}
void Prepare() {
g[1] = g[2] = 1;
f[1] = f[2] = 1;
invf[1] = invf[2] = 1;
for (int i = 3; i <= kN; ++ i) {
g[i] = 1;
f[i] = (f[i - 1] + f[i - 2]) % mod;
invf[i] = QPow(f[i], mod - 2);
}
Euler();
for (int d = 1; d <= kN; ++ d) {
for (int j = 1; d * j <= kN; ++ j) {
if (mu[j] == 1) {
g[d * j] = g[d * j] * f[d] % mod;
} else if (mu[j] == -1) {
g[d * j] = g[d * j] * invf[d] % mod;
}
}
}
invp[0] = prod[0] = 1;
for (int i = 1; i <= kN; ++ i) {
prod[i] = prod[i - 1] * g[i] % mod;
invp[i] = QPow(prod[i], mod - 2);
}
}
//=============================================================
int main() {
Prepare();
int T = read();
while (T -- ){
n = read(), m = read(), ans = 1;
if (n < m) std::swap(n, m);
for (LL l = 1, r = 1; l <= m; l = r + 1) {
r = std::min(n / (n / l), m / (m / l));
ans = (ans * QPow(prod[r] * invp[l - 1] % mod, 1ll * (n / l) * (m / l))) % mod;
}
printf("%lld\n", ans);
}
return 0;
}
P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题
给定参数 \(p\),有 \(T\) 组数据,每次给定参数 \(A,B,C\),求:
\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,k)}\right)^{f(type)} \]其中 \(f(type)\) 的取值如下:
\[f(type) = \begin{cases} 1 &type = 0\\ i\times j\times k &type = 1\\ \gcd(i,j,k) &type = 2 \end{cases}\]\(1\le A,B,C\le 10^5\),\(10^7\le p\le 1.05\times 10^9\),\(p\in \mathbb{P}\),\(T=70\)。
2.5S,128MB。
先化下原式,原式等于:
发现每一项仅与两个变量有关,设:
发现 \(\prod\) 可以随意交换,则原式等价于:
考虑在 \(type\) 取值不同时,如何快速求得 \(f_1\) 与 \(f_2\)。
一共有 \(6\) 个需要推导的式子,不妨就叫它们 \(1\sim 6\) 面了(
type = 0
对于 1 面,显然有:
预处理阶乘 + 快速幂即可,单次计算时间复杂度 \(O(\log n)\)。
再考虑 2 面,套路地枚举 \(\gcd\),显然有:
指数是个套路,可以看这里:P3455 [POI2007]ZAP-Queries。于是有:
代回原式,略做处理,则原式等于:
像「SDOI2017」数字表格 一样,考虑枚举 \(t=kd\) 和 \(d\),则原式等于:
设:
线性筛预处理 \(\mu\) 后,\(g_0(t)\) 可以用埃氏筛预处理,时间复杂度 \(O(n\log n)\)。再代回原式,原式等于:
预处理 \(g_0(t)\) 的前缀积和前缀积的逆元,复杂度 \(O(n\log n)\)。
数论分块 + 快速幂计算即可,单次时间复杂度 \(O(\sqrt n\log n)\)。
type = 1
考虑 3 面,把 \(\prod k\) 扔到指数位置,有:
再把 \(\prod j\) 也扔到指数位置,引入 \(\operatorname{sum}(n) = \sum_{i=1}^{n} i = \frac{n(n+1)}{2}\),原式等于:
预处理 \(i^i\) 的前缀积,复杂度 \(O(n\log n)\)。
指数可以 \(O(1)\) 算出,然后快速幂,单次时间复杂度 \(O(\log n)\)。
根据费马小定理,指数需要对 \(p - 1\) 取模。注意 \(p-1\) 不是质数,计算 \(\operatorname{sum}\) 时不能用逆元,但乘不爆 LL
,直接算就行。
再考虑 4 面,发现 \(k\) 与 \(\gcd\) 无关,则同样把 \(\prod k\) 扔到指数位置,则有:
套路地枚举 \(\gcd\),原式等于:
大力化简指数,有:
指数化不动了,代回原式,原式等于:
同 2 面的情况,先展开一下,再枚举 \(t=kd\) 和 \(d\),原式等于:
与二面相同,设:
\(g_1(t)\) 可以用埃氏筛套快速幂筛出,时间复杂度 \(O(n\log^2 n)\)。再代回原式,原式等于:
同样预处理 \(g_1(t)\) 的前缀积及其逆元,时间复杂度 \(O(n\log n)\)。
整除分块 + 快速幂即可,单次时间复杂度 \(O(\sqrt n\log n)\)。
注意指数的取模。
type = 2
考虑 5 面,手段同上,大力反演化简一波,再调换枚举对象,则有:
引入 \(\operatorname{fac}(n) = \prod_{i=1}^{n} i\),再根据枚举对象调整一下指数,原式等于:
指数中出现了一个经典的狄利克雷卷积的形式,对其进行反演。
将 \((\operatorname{Id}\ast \mu) (n)= \varphi (n)\) 代入原式,原式等于:
预处理 \(t^{\varphi(t)}\) 的前缀积及逆元,阶乘的前缀积及阶乘逆元,\(\pmod {p-1}\) 下的 \(\varphi(t)\) 的前缀和(指数
),时间复杂度 \(O(n\log n)\)。
同样整除分块 + 快速幂即可,单次时间复杂度 \(O(\sqrt n\log n)\)。
然后是最掉 sans 的 6 面。有 \(\gcd(i,j,k) = \gcd(\gcd(i,j), k)\),考虑先枚举 \(\gcd(i,j)\),然后套路化式子,则有:
先考虑最外面的指数,这也是个套路,可以参考 一个例子。用 \(\operatorname{Id} = \varphi \ast 1\) 反演,显然有:
再考虑里面的指数,发现这式子在 2 面已经推了一遍了,于是直接拿过来用,有:
将化简后的两个指数代入原式,原式等于:
与 2、4 面同样套路地,考虑枚举 \(t=yd\) 和 \(d\),再略作调整,原式等于:
发现要同时枚举 \(d\) 和 \(x\),化不动了。
从题解里学到一个比较神的技巧,考虑把 \(d\) 拆成 \(x\) 和 \(\frac{d}{x}\) 分别计算贡献再相乘,即分别计算下两式:
先考虑 \(x\) 的情况,首先把枚举 \(x\) 调整到最外层。设 \(\operatorname{lim}=\max(a,b,c)\),则原式等于:
把 \(\prod {t}\) 挪到指数位置,原式等于:
指数中又出现了一个经典的狄利克雷卷积的形式,对其进行反演。
将 \((\mu \ast 1) (n)= \epsilon (n)=[n=1]\) 代入原式,原式等于:
得到了一个非常优美的式子,而且发现这个式子是 5 面最终答案的一部分。同 5 面的做法,直接整除分块即可。
再考虑 \(\frac{d}{x}\) 的情况,同上先把枚举 \(x\) 放到最外层,并调整一下指数,则原式等于:
考虑枚举 \(dx\),替换原来的 \(d\),注意一下这里的倍数关系。原式等于:
发现最内层的式子 \(\prod_{d|t}d^{\mu\left(\frac{t}{d}\right)}\),即为二面处理过的 \(g_0(t)\)。直接代入,原式等于:
一个小结论,证明可以看 这里:
则原式等于:
于是可以先对外层整除分块,再对内层整除分块。
然后就做完了,哈哈哈。
一些实现上的小技巧:
- 逆元能预处理就预处理。
- 注意对指数取模,模数为 \(p-1\)。
//知识点:莫比乌斯反演
/*
By:Luckyblock
用了比较清晰易懂的变量名,阅读体验应该会比较好。
vsc 的自动补全真是太好啦!
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
using std::min;
using std::max;
#define LL long long
const int Lim = 1e5;
const int kN = 1e5 + 10;
//=============================================================
LL A, B, C, mod, ans;
int T, p_num, p[kN];
bool vis[kN];
LL mu[kN], phi[kN], fac[kN], g[2][kN];
LL sumphi[kN], prodt_phi[kN], prodi_i[kN], prodg[2][kN];
LL inv[kN], inv_fac[kN], inv_prodt_phi[kN], inv_prodg[2][kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) {
w = (w << 3) + (w << 1) + (ch ^ '0');
}
return f * w;
}
void Chkmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
x_ %= mod;
y_ %= mod - 1;
LL ret = 1;
for (; y_; y_ >>= 1ll) {
if (y_ & 1) ret = ret * x_ % mod;
x_ = x_ * x_ % mod;
}
return ret;
}
LL Inv(LL x_) {
return QPow(x_, mod - 2);
}
LL Sum(LL n_) {
return ((n_ * (n_ + 1ll)) / 2ll) % (mod - 1);
}
void Euler() {
vis[1] = true, mu[1] = phi[1] = 1; //初值
for (int i = 2; i <= Lim; ++ i) {
if (! vis[i]) {
p[++ p_num] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for (int j = 1; j <= p_num && i * p[j] <= Lim; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
phi[i * p[j]] = phi[i] * p[j];
break;
}
mu[i * p[j]] = -mu[i];
phi[i * p[j]] = phi[i] * (p[j] - 1);
}
}
}
void Prepare() {
Euler();
inv[1] = fac[0] = prodt_phi[0] = prodi_i[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
g[0][i] = g[1][i] = 1;
fac[i] = 1ll * fac[i - 1] * i % mod;
sumphi[i] = (sumphi[i - 1] + phi[i]) % (mod - 1);
prodi_i[i] = prodi_i[i - 1] * QPow(i, i) % mod;
if (i > 1) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
prodt_phi[i] = prodt_phi[i - 1] * QPow(i, phi[i]) % mod;
inv_prodt_phi[i] = Inv(prodt_phi[i]);
}
for (int d = 1; d <= Lim; ++ d) {
for (int j = 1; d * j <= Lim; ++ j) {
int t = d * j;
if (mu[j] == 1) {
g[0][t] = g[0][t] * d % mod;
g[1][t] = g[1][t] * QPow(1ll * d, 1ll * t * t) % mod;
} else if (mu[j] == -1) {
g[0][t] = g[0][t] * inv[d] % mod;
g[1][t] = g[1][t] * Inv(QPow(1ll * d, 1ll * t * t)) % mod;
}
}
}
inv_prodg[0][0] = prodg[0][0] = 1;
inv_prodg[1][0] = prodg[1][0] = 1;
inv_prodt_phi[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
for (int j = 0; j <= 1; ++ j) {
prodg[j][i] = prodg[j][i - 1] * g[j][i] % mod;
inv_prodg[j][i] = Inv(prodg[j][i]);
}
}
}
LL f1(LL a_, LL b_, LL c_, int type) {
if (! type) return QPow(fac[a_], b_ * c_);
if (type == 1) return QPow(prodi_i[a_], Sum(b_) * Sum(c_));
LL ret = 1, lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(fac[a_ / l], (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
return ret;
}
LL f2_2(LL a_, LL b_) {
LL ret = 1;
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
ret = ret * QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l)) % mod;
}
return ret;
}
LL f2(LL a_, LL b_, LL c_, int type) {
LL ret = 1;
if (! type) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l));
ret = (ret * QPow(val, c_)) % mod;
}
} else if (type == 1) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[1][r] * inv_prodg[1][l - 1], Sum(a_ / l) * Sum(b_ / l));
ret = (ret * QPow(val, Sum(c_))) % mod;
}
} else {
LL lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(f2_2(a_ / l, b_ / l), (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
}
return ret;
}
//=============================================================
int main() {
T = read(), mod = read();
Prepare();
while (T -- ) {
A = read(), B = read(), C = read();
for (int i = 0; i <= 2; ++ i) {
ans = f1(A, B, C, i) * f1(B, A, C, i) % mod;
ans = ans * Inv(f2(A, B, C, i)) % mod * Inv(f2(A, C, B, i)) % mod;
printf("%lld ", ans);
}
printf("\n");
}
return 0;
}
写在最后
参考资料:
Oi-Wiki-莫比乌斯反演
算法学习笔记(35): 狄利克雷卷积 By: Pecco
题解 SP5971 【LCMSUM - LCM Sum】 - BJpers2 的博客
题解 SP5971 【LCMSUM - LCM Sum】 - Venus 的博客
题解 P3327 【[SDOI2015]约数个数和】 - suncongbo 的博客
把 Oi-Wiki 上的内容进行了复制 整理扩充。
我是个没有感情的复读机(大雾)