莫比乌斯函数入门
缘起
入坑数论中的重要的部分——莫比乌斯函数
分析
本文的例题来自 Hdu 6053 TrickGCD 莫比乌斯函数的入门题.
给定一个长度为n的数组A， zhu同学想知道有多少个不同的数组 B 满足以下条件
1. 1 <= B_i <= A_i
2. 对任何满足 1<=l<=r<=n 的(l, r), 我们有 gcd(b_l,...,b_r)>=2
【输入】
多样例, 第一行是样例数, 每个样例开始于 n, 然后一行是包含 n 个整数Ai, 1<=n, Ai<=1e5
【输出】
对第 k 个询问, 输出 Case #k , 然后输出答案, 因为答案太大了, 所以答案要 mod掉 1e9+7
【样例输入】
1
4
4 4 4 4
【样例输出】
Case #1: 17
【限制】
2500ms 256MB
【来源】
2017暑期多校训练第二场
要解决这道题目，它要求莫比乌斯（Mobius）函数作为知识背景. 所以我们先来学习一下mobius函数.
mobius函数是一个数论函数，它是由德国著名天文学家以及数学家莫比乌斯（没错，就是发现莫比乌斯环这种反人类的拓扑结构的天才人物），它的定义如下
mobius函数是一个积性函数, 即
这里顺便说一下，如果把 gcd(m,n)=1这个条件去掉的话，则就是完全积性函数的定义.
既然介绍了mobius函数的定义，那么我们自然要琢磨如何计算mobius函数. 但是这种计算分成两类
- 给定 m, 要你计算
- 给定 m, 要你一口气计算 [1,...,m] 中所有整数的mobius函数值
其实这种事情我们在欧拉函数（也是一个数论函数）的时候也做过. 这里对于 mobius函数也给出这两类问题的模板.
问题1
显然根据 mobius 函数的定义, 我们只需要将 m 分解素因子然后按mobius函数的定义直接撸就可以了.
ili mobius(int m)
{
if (m == 1) return 1;
int tot = 0; // tot 是 m 不同素因子的个数
for (re i = 2, num; i * i <= m; i++)
{
if (m % i) continue;
++tot;
num = 0; // num 是 m 的 i 素因子的个数
while (m % i == 0)
{
m /= i;
if (++num > 1) return 0;
}
}
return tot & 1 ? -1 : 1;
}
问题2
在线性筛的过程中顺便计算mobius函数即可.
先上模板再解释该模板.
上面过程本质上来讲就是
假设已经遍历到了
（
）, 然后已经发现了素数 prime[1,..,tot]，则就用
去筛掉一些合数（即代码的第16行）并且计算相应的 mobius 值.
那么我们至少会有两个疑问
疑问1: 23行为什么当mobius函数值为0的话就结束本轮的for循环? 因为原本你可以使用 {i * prime[1], i * prime[2],...i * prime[tot]} 来筛掉更多的合数呀~ 但是现在你只使用了 {i * prime[1], i * prime[2],..., i * prime[j]}, 1 <= j <= tot 来筛掉合数啊~
疑问2: 其实是由疑问1产生的，因为你现在使用的并不是完整的 {i * prime[1], i * prime[2],...i * prime[tot]} 来筛掉合数. 所以又凭什么断定第8行 !isp[i] 成的话 i 就一定是素数呢?
鉴于疑问2由疑问1产生，所以这两个疑问我们一起来回答.
或者更确切的说，我们试图说明上面的mobius算法是正确的.
我们称
是由
发现的. 如果上面的伪代码的 第15行发生了.
做一般性的考虑, 令
，其中, {
} 是两两不同的质数且
那么根据伪代码第24行，
一定是由
发现的. 所以我们就知道上述伪代码的每一次内层for循环发现的合数集合的交集都是空集. 这是因为对于任何
,
的
是唯一的. 更确切讲,
不可能由
发现. 因为第13行是按照 prime[1]、prime[2]、...prime[tot] 这个从小到大的顺序遍历已经找到的素数的. 在遍历到素数
的时候属于
的内层for循环就已经 break掉了（因为
）.
这是理解上述模板的核心观察.
好, 现在基于该核心观察来回答疑问 1 和 疑问 2. 首先，任何一个[1,..,m]中的合数都会被筛掉，或者说，任何一个合数 都会被上述过程发现. 因为 合数, 可以做算术基本分解 , 那么根据核心观察， 将被 i = p_1^{r_1-1}p_2^{r_2}...*p_k^{r_k} 发现所以一定会在第行的循环（即外层循环）为i的时候把t 给发现出来. 其次，每个合数都只会恰好被发现一次. 这一点上面已经说了, 不再赘述. 然后我们回答为什么第8行发生 ! isp[i] 的时候, 一定是素数? 反证法，如果 是合数，令 , {} 是严格小于 的质数, 则根据核心观察（或者说就是上述伪代码的算法），知道 一定尚未被被发现. 进而 一定尚未被发现,... 然后不断的减少指数幂, 最后我们会推出 是尚未识别的质数. 这是不可能的. 因为
于是，我们讲清楚了mobius函数的线性筛算法. 该算法的复杂度是 O(n) 的.
有了mobius函数的基础知识，现在我们来看看本题. 题目要求的其实是计数
那么，自然的想法就是枚举
的值. 这种枚举公因子的手法在容斥原理的题目中多次使用过了.
而因为
，所以 d 在 [2, 1e5] 范围内 . d 对答案的贡献简单的使用乘法原理便可求出
特别的，如果
的话, gx(d) = 0. 根据容斥原理
那么我们就要求出 d至少是 i 个素数乘积的情况 这一计数结果. 而这是可以通过问题2的模板 O(n) 预处理出来的. 即通过mobius函数的线性筛法我们可以知道 [1,1e5]内所有的数的 mobius 函数值. 对于是 i 个不同素数乘积的数，它的mobius值恰好就是
, 而这刚好就是容斥原理的每一项的系数的相反数.
于是，我们可以切一版代码出来
//#include "stdafx.h"
//#define LOCAL
#pragma GCC optimize(2)
#pragma G++ optimize(2)
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <iostream>
#include <string>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
using namespace std;
#define int long long
#define re register int
#define ci const int
typedef pair<int, int> P;
#define FE(cur) for(re h = head[cur], to, len; ~h; h = g[h].nxt)
#define ilv inline void
#define ili inline int
#define ilc inline char
#define ild inline double
#define ilp inline P
#define LEN(cur) (hjt[cur].r - hjt[cur].l)
#define MID(cur) (hjt[cur].l + hjt[cur].r >> 1)
typedef vector<int>::iterator vit;
typedef set<int>::iterator sit;
typedef map<int, int>::iterator mit;
namespace fastio
{
const int BUF = 1 << 21;
char fr[BUF], fw[BUF], *pr1 = fr, *pr2 = fr;int pw;
ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) + fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1++; }
ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
ilv pc(char c) { if (pw >= BUF) flush(); fw[pw++] = c; }
ili read(int &x)
{
x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
while(!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
while(isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = gc();
x *= f; return 1;
}
ili read(double &x)
{
int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
while (isdigit(c)) { xx = (xx << 3) + (xx << 1) + (c ^ 48), c = gc(); }
x = xx;
if (c ^ '.') { x = f * xx; return 1; }
c = gc();
while (isdigit(c)) x += (c ^ 48) * (fraction /= 10), c = gc();
x *= f; return 1;
}
ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 + 48); }
ilv writeln(int x) { write(x);pc(10); }
ili read(char *x)
{
char c = gc(); if (!~c) return EOF;
while(!isalpha(c) && !isdigit(c)) c = gc();
while (isalpha(c) || isdigit(c)) *x++ = c, c = gc();
*x = 0; return 1;
}
ili readln(char *x)
{
char c = gc(); if (!~c) return EOF;
while(c == 10) c = gc();
while(c >= 32 && c <= 126) *x++ = c, c = gc();
*x = 0; return 1;
}
ilv write(char *x) { while(*x) pc(*x++); }
ilv writeln(char *x) { write(x); pc(10); }
ilv write(char c) { pc(c); }
ilv writeln(char c) { write(c); pc(10); }
} using namespace fastio;
const int maxn = 1e5+5, mod = 1e9+7, inf = 1ll<<60;
int n, a[maxn], mu[maxn], prime[maxn], isp[maxn], tot, cmin;
ilv mobius(int m)
{
mu[1] = 1;
for (re i = 2; i<= m; i++)
{
if (!isp[i])
{
prime[++tot] = i;
mu[i] = -1;
}
for (re j = 1, t; j <= tot && i * prime[j] <= m; j++)
{
t = prime[j] * i;
isp[t] = 1;
if (i % prime[j])
{
mu[t] = -mu[i];
}
else
{
mu[t] = 0;
break;
}
}
}
}
ili gx(int d)
{
int ans = 1;
for (re i = 1; i <= n; i++)
{
ans *= a[i] / d;
ans %= mod;
}
return ans;
}
ili kk()
{
int ans = 0;
for (re d = 2; d <= cmin; d++)
{
if (mu[d])
{
ans -= mu[d] * gx(d); // -mu[d] 才是容斥原理中的系数
ans = (ans + mod) % mod;
}
}
return ans;
}
signed main()
{
#ifdef LOCAL
freopen("d:\\data.in", "r", stdin);
// freopen("d:\\my.out", "w", stdout);
#endif
mobius(100000);
int kase; read(kase);
for (re i = 1; i <= kase; i++)
{
read(n);
cmin = inf;
for (re j = 1; j <= n; j++) read(a[j]), cmin = min(cmin, a[j]);
write("Case #"), write(i), write(": "), writeln(kk());
}
flush();
return 0;
}
但是遗憾的是，上面的代码被T掉了. 被T掉的原因也是很显然的，因为 gx 的复杂度是 O(n) 啊, 所以跑一个case的最坏复杂度达到了
，
那么怎么优化呢? 也就是要优化
的计算. 不要O(n) 计算. 咱们换个角度来看这个问题，其实
不就是一堆的 1、2、3、.... 么? 所以
注意，上述形式里面已经出现了幂次, 也就是我们可以使用快速幂来优化. 在此之前，我们需要知道
、
、...是多少，而我们很容易知道
，而要得到数组在一个权值范围内的元素的个数只需要简单的哈希就行了
于是此题就破了~
//#include "stdafx.h"
//#define LOCAL
#pragma GCC optimize(2)
#pragma G++ optimize(2)
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <iostream>
#include <string>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
using namespace std;
#define int long long
#define re register int
#define ci const int
typedef pair<int, int> P;
#define FE(cur) for(re h = head[cur], to, len; ~h; h = g[h].nxt)
#define ilv inline void
#define ili inline int
#define ilc inline char
#define ild inline double
#define ilp inline P
#define LEN(cur) (hjt[cur].r - hjt[cur].l)
#define MID(cur) (hjt[cur].l + hjt[cur].r >> 1)
typedef vector<int>::iterator vit;
typedef set<int>::iterator sit;
typedef map<int, int>::iterator mit;
namespace fastio
{
const int BUF = 1 << 21;
char fr[BUF], fw[BUF], *pr1 = fr, *pr2 = fr;int pw;
ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) + fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1++; }
ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
ilv pc(char c) { if (pw >= BUF) flush(); fw[pw++] = c; }
ili read(int &x)
{
x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
while(!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
while(isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = gc();
x *= f; return 1;
}
ili read(double &x)
{
int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
while (isdigit(c)) { xx = (xx << 3) + (xx << 1) + (c ^ 48), c = gc(); }
x = xx;
if (c ^ '.') { x = f * xx; return 1; }
c = gc();
while (isdigit(c)) x += (c ^ 48) * (fraction /= 10), c = gc();
x *= f; return 1;
}
ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 + 48); }
ilv writeln(int x) { write(x);pc(10); }
ili read(char *x)
{
char c = gc(); if (!~c) return EOF;
while(!isalpha(c) && !isdigit(c)) c = gc();
while (isalpha(c) || isdigit(c)) *x++ = c, c = gc();
*x = 0; return 1;
}
ili readln(char *x)
{
char c = gc(); if (!~c) return EOF;
while(c == 10) c = gc();
while(c >= 32 && c <= 126) *x++ = c, c = gc();
*x = 0; return 1;
}
ilv write(char *x) { while(*x) pc(*x++); }
ilv writeln(char *x) { write(x); pc(10); }
ilv write(char c) { pc(c); }
ilv writeln(char c) { write(c); pc(10); }
} using namespace fastio;
const int maxn = 1e5+5, mod = 1e9+7, inf = 1ll<<60;
int n, a[maxn], mu[maxn], prime[maxn], isp[maxn], tot, cmin, cmax, cnt[maxn << 1], pref[maxn << 1];
ilv mmax(int &a, int b)
{
if (a < b)
{
a = b;
}
}
ilv mmin(int &a, int b)
{
if (a > b)
{
a = b;
}
}
ilv mobius(int m)
{
mu[1] = 1;
for (re i = 2; i<= m; i++)
{
if (!isp[i])
{
prime[++tot] = i;
mu[i] = -1;
}
for (re j = 1, t; j <= tot && i * prime[j] <= m; j++)
{
t = prime[j] * i;
isp[t] = 1;
if (i % prime[j])
{
mu[t] = -mu[i];
}
else
{
mu[t] = 0;
break;
}
}
}
}
ili ksm(int a, int b)
{
int ans = 1;
while (b)
{
if (b & 1)
{
ans = ans * a;
ans %= mod;
}
if (b > 1)
{
a = a * a;
a %= mod;
}
b >>= 1;
}
return ans;
}
ili gx(int d)
{
int ans = 1;
for (re i = 1; i * d <= cmax; i++) // [i * d, (i + 1) * d) 范围
{
ans *= ksm(i, pref[(i + 1) * d - 1] - pref[i * d - 1]); // 快速幂优化gx的计算
ans %= mod;
}
return ans;
}
ili kk()
{
int ans = 0;
for (re d = 2; d <= cmin; d++)
{
if (mu[d])
{
ans -= mu[d] * gx(d);
ans = (ans + mod) % mod;
}
}
return ans;
}
signed main()
{
#ifdef LOCAL
freopen("d:\\data.in", "r", stdin);
freopen("d:\\my.out", "w", stdout);
#endif
mobius(1e5);
int kase; read(kase);
for (re i = 1; i <= kase; i++)
{
read(n);
cmax = 0, cmin = inf;
memset(cnt, 0, sizeof(cnt));
for (re j = 1; j <= n; j++) read(a[j]),mmin(cmin, a[j]), mmax(cmax, a[j]), ++cnt[a[j]];
for (re i = 1; i <= cmax << 1; i++) pref[i] = pref[i - 1] + cnt[i];
write("Case #"), write(i), write(": "), writeln(kk());
}
flush();
return 0;
}
ac情况
accepted 280ms 8988kB
小结
mobius函数的作用就在于它O(n)时间预处理出了所有恰好由若干两两不同的质数相乘得到的数. 而这些数恰好是我们在容斥原理中要用到的数.