入坑數論中的重要的部分——莫比烏斯函數
分析本文的例題來自 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, 要你一口氣計算 [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函數即可.
先上模板再解釋該模板.
上面過程本質上來講就是
假設已經遍歷到了
那麼我們至少會有兩個疑問
疑問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算法是正確的.
我們稱
做一般性的考慮, 令 ，其中, {
那麼根據偽代碼第24行，
這是理解上述模板的核心觀察.
好, 現在基於該核心觀察來回答疑問 1 和 疑問 2. 首先，任何一個[1,..,m]中的合數都會被篩掉，或者說，任何一個合數
於是，我們講清楚了mobius函數的線性篩算法. 該算法的複雜度是 O(n) 的.
有了mobius函數的基礎知識，現在我們來看看本題. 題目要求的其實是計數
那麼，自然的想法就是枚舉
而因為
特別的，如果
那麼我們就要求出 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) 計算. 咱們換個角度來看這個問題，其實
注意，上述形式裡面已經出現了冪次, 也就是我們可以使用快速冪來優化. 在此之前，我們需要知道
於是此題就破了~
//#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)時間預處理出了所有恰好由若干兩兩不同的質數相乘得到的數. 而這些數恰好是我們在容斥原理中要用到的數.