线性筛与莫比乌斯反演

目录

• 线性筛与莫比乌斯反演
  • 线性筛
    • 代码
    • 讲解
  • 莫比乌斯反演
    • 定理
    • 代码
    • 常见的定理
    • 莫比乌斯反演的证明
  • 一些例题（难题）
    • 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 mod 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} ]

莫比乌斯反演的证明

证明：

[egin{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) ]

题解

一个莫比乌斯反演然后化式子。

[egin{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} floor} sum _{j=1}^{lfloor frac{m}{d} floor} ij [gcd(i,j)=1]\ end{aligned} ]

这个就是一个更换枚举相的操作了，是个套路。你先枚举所有可能的(gcd)再计算这种(gcd)的贡献。

比如前面的那个(d)就是我们枚举的(gcd)，后面所有可能的数对，就是在(lfloor frac{n}{d} floor)和(lfloor frac{m}{d} floor)中的所有互质的数对的乘积在乘上(d)。

这个可以简单理解一下，就是两个数分别除以他们的最大公因数，然后两个数肯定是互质的。
但其对于答案的贡献就多除以了一个(d)，所以要乘回来。

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

这个就是运用了前面的公式(sum limits _{d|n} mu(d)=[n=1])来替代了([gcd(i,j)=1])的条件。（这个就是套路了）

然后我们继续推：

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

这个也是套路，把(x)提前了。就是改成了枚举(x)看看它的对于答案的贡献是多少。

很容易发现，就是在([1,lfloor frac{n}{dx} floor])中的所有数乘上(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))，(mod)的常数还很大）所以就需要来用套路的整除分块了。
就是把后面两个 (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