線性篩與莫比烏斯反演

時間 2019-11-09

原文原文鏈接

目錄c++

線性篩與莫比烏斯反演

線性篩與莫比烏斯反演

和上篇文章同樣，一直沒有研究這個東西，結果又考了GG……TAT
下定決心學一學，搞好這個東西。git

線性篩

篩質數有不少方法，好像很厲害的有洲閣篩 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 、杜教篩 \(O(n^{\frac{2}{3}})\)（~~然而我都不會QAQ~~）(其實這些不是用來篩素數的2333 用來篩積性函數前綴和的 \(Update \ 2018 \cdot 4 \cdot 9\) )，還有暴力篩（就是枚舉一個數的倍數）複雜度是 \(O(n \ln n)\) 的。
我只學了比較簡單並且實用的線性篩法。
這種篩法是避免一個數被重複篩幾遍，因此效率均攤下來能夠達到線性。(網上有證實)github

代碼

const int N = 100000;
bool is_prime[N+100];
int prime[N], cnt = 0;
void find_prime() {
    Set(is_prime, true);
    is_prime[0] = is_prime[1] = false;
    For(i, 2, N) {
        if (is_prime[i]) prime[++cnt] = i;
    For(j, 1, cnt) {
        if (i * prime[j] > N) break;
        is_prime[i * prime[j]] = false;
        if (i % prime[j] == 0) break; //here
    }
    }
}

講解

這個代碼有一個關鍵點就是上面的\(here\) 這個意義就是對於一個合數\(m\)能夠分解爲\(m=p_1^{r_1}*...*p_n^{r_n}\)其中
\(p_i\)爲質數，那麼咱們篩\(m\)的時候以前把\(p_1\)篩掉了，因此在枚舉\(i\)的時候。函數

若是\(i\)爲素數沒問題，直接向後繼續推（由於篩出的質數都相似\(m=p_1*p_2\)的形式，因此不可能重複）。
若是爲合數，那麼\(i\)能夠分解成\(i={p_1}^{r_1}*...*{p_n}^{r_n}\)形式其中\(p_1-p_n\)是遞增的，
那麼\(p_1\)是最小的那個質數。\(i \bmod p_1 = 0\)的時候，就不用繼續枚舉了，因此咱們就只能篩出不大於\(p_1\)的質數\(*i\)。

複雜度好像是線性的（ \(O(n)\) ），我不太會證複雜度。。優化

莫比烏斯反演

莫比烏斯反演不少時候都能大大簡化運算……spa

定理

\(F(n)\)和\(f(n)\)是定義在非負整數集合上的兩個函數，而且知足條件\[F(n)=\sum \limits \limits _{d|n}{f(d)}\]。那麼咱們就能獲得結論：
\[ f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d}) \]code

在上面的公式中有一個\(\mu(d)\)函數（莫比烏斯函數），它的定義以下：ci

若\(d=1\),那麼\(\mu(d)=1\)
若\(d=p_1p_2...p_k\)，\(p_i\)均爲互異質數，那麼\(\mu(d)=(-1)^{k}\)。這個個人理解就是\(d\)的質因數個數爲偶數的話，那麼\(\mu(d)=1\) 不然爲 \(-1\) 。get
其餘狀況下\(\mu(d)=0\)這個就是對上面那條的拓展了，就是指的\(d\)沒有一個平方因子，或者說沒有一個質因子的次數大於 \(1\) 。博客

代碼

const int N = 100100;
bool is_prime[N+100];
int mu[N+100] = {0, 1}, cnt = 0, prime[N+100];
void init() {
    Set(is_prime, true);
    is_prime[1] = false;
    For (i, 2, N) {
        if (is_prime[i]) {
            prime[++cnt] = i;
            mu[i] = -1; //質數的質因子個數確定爲奇數個就是1
        }
        For (j, 1, cnt) {
            if (i * prime[j] > N) break;
            is_prime[i * prime[j]] = false;
            if (i % prime[j]) mu[i * prime[j]] = -mu[i]; //多了一個質因子直接變爲原來結果的相反數
            else {
                mu[i * prime[j]] = 0; //這個將要被篩的數至少具備兩個prime[j]的因子
                break;
            }
        }
    }
}

常見的定理

而後還要提一下的就是一些常見的定理，證實嘛……

通常都是先分解質因數，而後再根據組合數性質去算，好比第一個。

要麼就是對於一些常見的反演格式進行反演，好比第二個。

\[ \sum _{d|n} \mu(d)=[n=1] \]

\[ \sum _{d|n} \frac{\mu(d)}{d}=\frac{\varphi(n)}{n} \]

莫比烏斯反演的證實

證實：

\[ \begin{aligned} \sum_{d|n}\mu(d)F(\frac{n}{d})&=\sum_{d|n}\mu(d)\sum \limits_{d'|\frac{n}{d}}f(d')\\ &=\sum_{d'|n}f(d')\sum_{d|\frac{n}{d'}}\mu(d)\\&=f(n) \end{aligned} \]

Q.E.D

一些例題（難題）

Luogu 【P1829】[國家集訓隊]Crash的數字表格

題意

\[ \sum_{i=1}^{n} \sum_{i=1}^{m} \mathrm{lcm}(i,j) \ (n,m \le 10^7) \]

題解

一個莫比烏斯反演而後化式子。

\[ \begin{aligned} ans &= \sum _{i=1}^{n} \sum _{i=1}^{m} \mathrm{lcm}(i,j)\\ &= \sum _{i=1}^{n} \sum _{i=1}^{m} \frac{ij}{\gcd(i,j)}\\ &= \sum _{d=1}^{min(n,m)} d \sum _{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum _{j=1}^{\lfloor \frac{m}{d} \rfloor} ij \ [\gcd(i,j)=1]\\ \end{aligned} \]

這個就是一個更換枚舉相的操做了，是個套路。你先枚舉全部可能的\(\gcd\)再計算這種\(\gcd\)的貢獻。

好比前面的那個\(d\)就是咱們枚舉的\(\gcd\)，後面全部可能的數對，就是在\(\lfloor \frac{n}{d} \rfloor\)和\(\lfloor \frac{m}{d} \rfloor\)中的全部互質的數對的乘積在乘上\(d\)。

這個能夠簡單理解一下，就是兩個數分別除以他們的最大公因數，而後兩個數確定是互質的。
但其對於答案的貢獻就多除以了一個\(d\)，因此要乘回來。

\[ \begin{aligned} ans =\sum_ {d=1}^{\min(n,m)} d \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{d} \rfloor} \sum_{x|\gcd(i,j)}\mu(x) * i * j \end{aligned} \]

這個就是運用了前面的公式\(\sum \limits _{d|n} \mu(d)=[n=1]\)來替代了\([gcd(i,j)=1]\)的條件。（這個就是套路了）

而後咱們繼續推：

\[ \begin{aligned} ans &=\sum_{d=1}^{\min(n,m)} d \sum_{x=1}^{min(\lfloor \frac{n}{d} \rfloor,\lfloor \frac{m}{d} \rfloor)}\mu(x) \ x^2 \sum_{i=1}^{\lfloor \frac{n}{dx} \rfloor} i \sum_{j=1}^{\lfloor \frac{m}{dx} \rfloor} j \end{aligned} \]

這個也是套路，把\(x\)提早了。就是改爲了枚舉\(x\)看看它的對於答案的貢獻是多少。

很容易發現，就是在\([1,\lfloor \frac{n}{dx} \rfloor]\)中的全部數乘上\(x\)

就是原來可行的\(i\)。而後咱們就能夠根據這個來優化了。

前面那兩個\(\sum \limits\)就是\(O(n \ln \ n)\)（令\(n=\min(n,m)\)）的複雜度。（就是調和級數 \(H(n) = \sum_{i=1}^{n} \frac{n}{i}\)）。後面的那兩個，直接用等差數列求和公式\(O(1)\)算。

但這個仍然過不去...（\(O(1.66*10^8)\)，\(\bmod\)的常數還很大）因此就須要來用套路的整除分塊了。
就是把後面兩個 \(\sum \limits\) 不少同樣答案的地方一塊兒處理掉，因此對於那個 \(\mu(x) \ x^2\) 還要記一個前綴和。

總複雜度\(O(\sum \limits _{i=1}^{n} \sqrt{\frac{n}{i}})=O(pass)\)這個我也不會算。。會算的大佬私聊啊23333

代碼

#include <bits/stdc++.h>
#define For(i, l, r) for (int i = (l), _end_ = (int)(r); i <= _end_; ++i)
#define Fordown(i, r, l) for (int i = (r), _end_ = (int)(l); i >= _end_; --i)
#define Set(a, v) memset(a, v, sizeof(a))
using namespace std;

bool chkmin(int &a, int b) { return b < a ? a = b, 1 : 0; }
bool chkmax(int &a, int b) { return b > a ? a = b, 1 : 0; }

void File() {
#ifdef zjp_shadow
    freopen("P2154.in", "r", stdin);
    freopen("P2154.out", "w", stdout);
#endif
}

typedef long long ll;
ll n, m;
const int N = 1e7 + 1e3;
const ll Mod = 20101009, inv2 = (Mod + 1) / 2;
bool is_prime[N];
int mu[N], prime[N], cnt;
ll sum[N];

inline void add(ll &x, ll y) {
    x = ((x + y) % Mod + Mod) % Mod;
}

void Get_Mu(int maxn) {
    int res;
    Set(is_prime, true);
    is_prime[0] = is_prime[1] = false;
    mu[1] = 1;
    For(i, 2, maxn) {
        if (is_prime[i])
            prime[++cnt] = i, mu[i] = -1;
        For(j, 1, cnt) {
            if ((res = i * prime[j]) > maxn) break;
            is_prime[res] = false;
            if (i % prime[j]) mu[res] = -mu[i];
            else { mu[res] = 0; break; }
        }
    }
    For(i, 1, maxn) add(sum[i], sum[i - 1] + (ll)mu[i] * i * i % Mod);
}

inline ll fsum(ll x) { return x * (1 + x) % Mod * inv2 % Mod; }

int main() {
    File();
    cin >> n >> m;
    Get_Mu(min(n, m));
    ll ans = 0;
    For(d, 1, min(n, m)) {
        int n_ = n / d, m_ = m / d;
        For(x, 1, min(n_, m_)) {
            if (!mu[x]) continue;
            int Next = min(n_ / (n_ / x), m_ / (m_ / x));
            add(ans, 1ll * d * (sum[Next] - sum[x - 1]) % Mod * fsum(n_ / x) % Mod * fsum(m_ / x) % Mod);
            x = Next;
        }
    }
    cout << ans << endl;
    return 0;
}