线性筛与莫比乌斯反演

目录

• 线性筛与莫比乌斯反演
  • 线性筛
    • 代码
    • 讲解
  • 莫比乌斯反演
    • 定理
    • 代码
    • 常见的定理
    • 莫比乌斯反演的证实
  • 一些例题（难题）
    • Luogu 【P1829】[国家集训队]Crash的数字表格
    • BZOJ3994 [SDOI2015]约数个数和

线性筛与莫比乌斯反演

和上篇文章同样，一直没有研究这个东西，结果又考了GG……TAT
下定决心学一学，搞好这个东西。

线性筛

筛质数有不少方法，好像很厉害的有洲阁筛 \(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)\) 的。
我只学了比较简单并且实用的线性筛法。
这种筛法是避免一个数被重复筛几遍，因此效率均摊下来能够达到线性。(网上有证实)

代码

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\)的时候。

1. 若是\(i\)为素数没问题，直接向后继续推（由于筛出的质数都相似\(m=p_1*p_2\)的形式，因此不可能重复）。
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)\) ），我不太会证复杂度。。

莫比乌斯反演

莫比乌斯反演不少时候都能大大简化运算……

定理

\(F(n)\)和\(f(n)\)是定义在非负整数集合上的两个函数，而且知足条件\[F(n)=\sum \limits \limits _{d|n}{f(d)}\]。那么咱们就能获得结论：
\[ f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d}) \]

在上面的公式中有一个\(\mu(d)\)函数（莫比乌斯函数），它的定义以下：

1. 若\(d=1\),那么\(\mu(d)=1\)

2. 若\(d=p_1p_2...p_k\)，\(p_i\)均为互异质数，那么\(\mu(d)=(-1)^{k}\)。这个个人理解就是\(d\)的质因数个数为偶数的话，那么\(\mu(d)=1\) 不然为 \(-1\) 。

3. 其余状况下\(\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;
}

---

BZOJ3994 [SDOI2015]约数个数和

题意

求\[\sum \limits _{i=1}^{n} \sum \limits _{j=1}^{m} d(ij)\]，\(d(x)\)表示\(x\)的约数个数。(\(n,m \le 10^5\))

题解

一个反演。当初数学一本通没看懂，真是本垃圾书。

肖大佬博客一看就懂了2333。只是中间有一步化\([gcd(i,j)=1]\)仍是习惯变两步，容易理解些QwQ