博客食用更香

一、素数

1、欧拉筛

时间复杂度：$O(n)$

constexpr int N = 1E6;

std::vector<int> primes;
std::vector<bool> st;

void init(int n) {
    st.assign(n + 1, false);
    primes.clear();
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes.push_back(i);
        }
        for (auto p : primes) {
            if (p * i > n) {
                break;
            }
            st[i * p] = true;
            if (i % p == 0) break;
        }
    }
}

2、MR-素数测试

时间复杂度：$O(k\log n)$

i64 mul(i64 a, i64 b, i64 m) {
    return static_cast<__int128>(a) * b % m;
}
i64 power(i64 a, i64 b, i64 m) {
    i64 res = 1 % m;
    for (; b; b >>= 1, a = mul(a, a, m))
        if (b & 1)
            res = mul(res, a, m);
    return res;
}
bool isprime(i64 n) {
    if (n < 2)
        return false;
    static constexpr int A[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
    int s = __builtin_ctzll(n - 1);
    i64 d = (n - 1) >> s;
    for (auto a : A) {
        if (a == n)
            return true;
        i64 x = power(a, d, n);
        if (x == 1 || x == n - 1)
            continue;
        bool ok = false;
        for (int i = 0; i < s - 1; ++i) {
            x = mul(x, x, n);
            if (x == n - 1) {
                ok = true;
                break;
            }
        }
        if (!ok)
            return false;
    }
    return true;
}

3、Pollard-Rho分解

分解质因数，返回一个分解后的质因数数组

时间复杂度：$O(n^{\tfrac{1}{4}}\log n)$

std::vector<i64> factorize(i64 n) {
    std::vector<i64> p;
    std::function<void(i64)> f = [&](i64 n) {
        if (n <= 10000) {
            for (int i = 2; i * i <= n; ++i)
                for (; n % i == 0; n /= i)
                    p.push_back(i);
            if (n > 1)
                p.push_back(n);
            return;
        }
        if (isprime(n)) {
            p.push_back(n);
            return;
        }
        auto g = [&](i64 x) {
            return (mul(x, x, n) + 1) % n;
        };
        i64 x0 = 2;
        while (true) {
            i64 x = x0;
            i64 y = x0;
            i64 d = 1;
            i64 power = 1, lam = 0;
            i64 v = 1;
            while (d == 1) {
                y = g(y);
                ++lam;
                v = mul(v, std::abs(x - y), n);
                if (lam % 127 == 0) {
                    d = std::gcd(v, n);
                    v = 1;
                }
                if (power == lam) {
                    x = y;
                    power *= 2;
                    lam = 0;
                    d = std::gcd(v, n);
                    v = 1;
                }
            }
            if (d != n) {
                f(d);
                f(n / d);
                return;
            }
            ++x0;
        }
    };
    f(n);
    std::sort(p.begin(), p.end());
    return p;
}

二、同余系

1、欧拉函数

1. 定义

对于正整数 $n$ ，欧拉函数是小于等于 $n$ 的正整数中与 $n$ 互质的数的数目，记作 $\varphi(n)$ 。特别的 $\varphi(1)=1$

即：$\varphi(n)=\sum\limits^n_{i=1}[\gcd(i,n)=1]$

在算数基本定理中，$N = p_1^{c_1} \times p_2^{c_2} \times ... \times p_k^{c_k}$，则：

$$\varphi(N)=N\times \frac{p_1-1}{p_1} \times \frac{p_2-1}{p_2}\times ...\times \frac{p_k-1}{p_k}=N\times \prod_{i=1}^k(1-\frac{1}{p}) $$

2. 求欧拉函数

首先欧拉函数是一个积性函数，当 $m,n$ 互质时，$\varphi(mn)=\varphi(m) \times \varphi(n)$

证明：把 $m$ 和 $n$ 分解质因数后带入欧拉函数的计算式

根据唯一分解定理： $n = p_1^{c_1} \times p_2^{c_2} \times ... \times p_k^{c_k}$

因此：$\varphi(n)=\varphi(p_1^{c_1}) \times ... \times \varphi(p_k^{c_k})$

对于任意一项 $\varphi(p_i^{c_i})=p_i^{c_i} - p_i^{c_i-1}$

证明：

从 $1$ 到 $p_i^{c_i}$ 中共有 $p_i^{c_i}$ 个数字

其中与 $p_i^{c_i}$ 不互质的有 $p_i,2p_i,3p_i,...,p_i\times p_i,...,p_i^{c_i-1}\times p_i$ 一共有 $p_i^{c-1}$ 项

根据容斥原理 $\varphi(p_i^{c_i})=p_i^{c_i} - p_i^{c_i-1}$

所以把公式再化简：$\varphi(p_i^{c_i})=p_i^{c_i} - p_i^{c_i-1}=p_i^{c_i}\times (1- \frac{1}{p_i})$

所以：

$$\begin{split} \varphi(n)=&\varphi(p_1^{c_1}) \times ... \times \varphi(p_k^{c_k})\\ =&(p_1^{c_1}-p_1^{c_1-1}) \times ... \times (p_k^{c_k} - p_k^{c_k-1}) \\ =&p_1^{c_1}\times (1- \frac{1}{p_1})\times p_2^{c_2}\times (1- \frac{1}{p_2}) \times ... \times p_k^{c_k}\times (1- \frac{1}{p_k})\\ =&p_1^{c_1}\times p_2^{c_2}\times...\times p_k^{c_k}\times(1- \frac{1}{p_1})\times (1- \frac{1}{p_2})\times...\times(1- \frac{1}{p_k}) \\ =&n\times \prod_{i=1}^k (1-\frac{1}{p_i}) \end{split} $$

3. 欧拉函数的性质

• 设 $p$ 为质数，若 $p|n$ 且 $p^2|n$ ，则 $\varphi(n)=\varphi(n/p)\times p$
• 设 $p$ 为质数，若 $p|n$ 且 $p^2\nmid n$ ，则 $\varphi(n)=\varphi(n/p)\times (p-1)$
• $\sum\limits_{d|n}\varphi(d)=n$

4. 代码实现

int phi(int n) {
    int res = n;
    for (int i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            while (n % i == 0) {
                n /= i;
            }
            res = res / i * (i - 1);
        }
    }
    if (n > 1) {
        res = res / n * (n - 1);
    }
    return res;
}

constexpr int N = 1E6;

std::vector<i64> phi, primes;
std::vector<bool> st;

void init(int n) {
    phi.resize(n);
    st.assign(n + 1, false);
    primes.clear();

	phi[1] = 1;
	for (int i = 2; i <= n; i++) {
		if (!st[i]) {
			primes.push_back(i);
			phi[i] = i - 1;
		}
		for (auto p : primes) {
			if (p * i > n) {
				break;
			}
			st[i * p] = true;
			if (i % p == 0) {
				phi[i * p] = phi[i] * p;
				break;
			}
			phi[i * p] = phi[i] * (p - 1);
		}
	}
}

2、扩展欧几里得算法

扩展欧几里得算法解决的问题：

1. 求解方程 $ax+by=\gcd(a,b)$ 的解

当 $b=0$ 时， $ax+by=a$ 故而 $x=1,y=0$

当 $b \ne 0$ 时

因为

$$\gcd(a,b)=\gcd(b,a\%b)$$

原式可写为：

$$bx^{\prime}+(a\%b)y^{\prime}=\gcd(b,a\%b)\\ bx^{\prime}+(a-\lfloor \frac ab \rfloor \times b)y^{\prime}=\gcd(b,a\%b)\\ ay^{\prime}+b(x^{\prime}-\lfloor \frac ab\rfloor \times y^{\prime})=\gcd(b,a\%b)=\gcd(a,b) $$

与原式比较得：

$$x=y^{\prime},y=x^{\prime}-\lfloor \frac ab \rfloor \times y^{\prime} $$

因此可以采取递归算法，先求出下一层的 $x^{\prime}$ 和 $y^{\prime}$ 再利用上述公式回代即可。

2. 求解更一般的方程 $ax+by=c$

前置知识：裴蜀定理

对于一个二元一次方程 $ax+by=c$，如果 $c$ 是 $gcd(a,b)$ 的倍数，那么这个方程一定有整数解，反之也成立。

证明：

充分性可通过exgcd求证，以证明过，下面证必要性，令 $k=gcd(a,b)$。

设 $a=a'k,b=b'k,c=c'k+r(r<k)$ 我们还不知道 $c$ 到底是不是 $k$ 的倍数，那就先设上个 $r$

原方程可变为：$a'kx+b'ky=c'k+r$

两边同时除以 $k$ ，得到 $a'x+b'y=c'+\frac{r}{k}$

看条件，有整数解，故 $\frac{r}{k}\in \Z$ 恒成立，$r$ 只能是 $0$ ，即 $c$ 一定是 $k$ 的倍数。

我们可以设 $d=\gcd(a,b)$ 则这个方程有解的条件是：当且仅当 $d \mid c$ 即 $d$ 整除 $c$ $(c \ \text{mod}\ d = 0)$

求解：

用扩展欧几里得求出 $ax_0+by_0=d$ 的解

则同时扩大 $\frac{c}{d}$ 倍： $a(x_0 \times \cfrac cd)+b(y_0 \times \cfrac cd)=c$

故而一个特解为：$x^{\prime}=x_0 \times \cfrac cd, \ y^{\prime}=y_0 \times \cfrac cd$

我们知道：非齐次方程的通解 = 非齐次方程的一个特解 + 齐次方程的通解

所以下一步求齐次方程：$ax+by=0$ 的通解

$$ax+by=0\\ ax=-by\\ x=-\frac{b}{a}y\\ 设d=\gcd(a,b)\\ \because a=a_1\times d\ ,b=b_1\times d\\ \therefore \gcd(a_1,b_1)=1\\ x=-\frac{b_1\times d}{a_1\times d}y\\ x=-\frac{b_1}{a_1}y\\ \because a_1,b_1互质，则x有整数解\\ a_1 \mid y\to y=k \times a_1 \to y=k\times \frac{a}{d}(k\in z)\\ \therefore 带入y得x=-k\frac{b}{d}\\ 同理y=x\times -\frac{a}{b}=k\frac{a}{d} $$

所以非齐次方程的通解为：

$$x=x^{\prime}-k\frac{b}{d},y=y^{\prime}+k\frac{a}{d},k\in z\\ $$

当然 $k\cfrac{b}{d}$ 前面的符号只需要保证在 $x$ 和 $y$ 的表达式中是相反的就可以

若令 $t=\cfrac bd$ ，则对于 $x$ 的最小非负整数解为 $(x^{\prime}\%t+t)\%t$

3. 求解一次同余方程 $ax\equiv b(\text{mod}\ m)$

等价于求：

$$ax=m\times (-y)+b\\ ax+my=b$$

有解条件为：$\gcd(a,m) \mid b$ ，然后用扩展欧几里得求解即可

特别的 当 $b=1$ 且 $a$ 与 $m$ 互质时，所求的 $x$ 即为 $a$ 的逆元

int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    retrun g;
}

3、逆元

乘法逆元的定义：若整数 $b,m$ 互质， 并且对于任意整数 $a$ ，如果满足 $b \mid a$ ，则存在一个 $x$ ，使得 $a \times b^{-1} \equiv a\times x(\text{mod}\ m)$ ，则称 $x$ 为 $b$ 的模 $m$ 的乘法逆元，记为 $b^{-1}(\text{mod}\ m)$。

1. 快速幂求逆元(m为质数)

推公式：

$$a\times b^{-1} \equiv a\times x(\text{mod} \ m)\\ 两边同时乘b得： a \equiv a \times b \times x(\text{mod} \ m)\\ 即 1\equiv b \times x(\text{mod} \ m)\\ 同 b \times x \equiv 1(\text{mod} \ m)\\ 由费马小定理可知，当m为质数时：\\ b^{m-1}\equiv 1(\text{mod} \ m)\\ 拆一个b出来\\ b\times b^{m -2}\equiv 1(\text{mod} \ m)\\ 对比上式得到 x=b^{m-2} $$

所以当 $m$ 为质数时，$b$ 的乘法逆元 $x=b^{m-2}(\text{mod} \ m)$ ，直接使用快速幂 qmi(b, m - 2, m)

2. 扩展欧几里得算法求逆元

前提：$a$ 有逆元的充要条件是 $a$ 与 $m$ 互质，所以有 $\gcd(a,m)=1$

假设 $a$ 的逆元为 $x$ ，那么有 $a\times x \equiv 1(\text{mod}\ m)$

等价：$ax+my=1$ 直接使用扩展欧几里得算法 exgcd(a, m, x, y)

如何保证求得的 $x$ 为正数？ cout<<(x + m) % m

3. 欧拉定理求逆元

欧拉定理：设 $m,a\in N^+$ ，且 $\gcd(a,m)=1$ ，则我们有 $a^{\varphi(m)}\equiv1(\text{mod}\ m)$ ，$\varphi(m)$ 为 $m$ 的欧拉函数

对上式定理进行变形得： $a \times a^{\varphi(m)-1}\equiv1(\text{mod}\ m)$

所以令 $a$ 的逆元是 $a^{-1}$ 即 $a^{-1}= a^{\varphi(m)-1}(\text{mod}\ m)$

4. 线性求逆元

有时候我们不一定是只求单个的逆元，需要求很多数的逆元，如果每次都用1、2、3方法复杂度就非常大了。

在线性时间复杂度求出 $[1,n]$ 在模 $p$ 意义下的逆元**（ $p$ 必须是质数）**

首先从 $p$ 入手， $p \div i=q\cdots r$ （其中 $q$ 是商， $r$ 是余数），移项得： $p=q\times i+r$

模上 $p$ 得：$0=q\times i+r(\text{mod}\ p)$

左右两边同时乘上 $i^{-1}\times r^{-1}$ 得： $0=q\times r^{-1} + i^{-1}(\text{mod}\ p)$

移项得：$i^{-1}=-q\times r^{-1}(\text{mod}\ p)$

又因为 $q=\lfloor \frac{p}{i} \rfloor$ ，$r=p\ \text{mod}\ i$ 带入上式得

$$i^{-1}=-\lfloor \frac{p}{i} \rfloor \times (p\ \text{mod}\ i)^{-1}(\text{mod}\ p)\\ i^{-1}=(-\lfloor \frac{p}{i} \rfloor \times (p\ \text{mod}\ i)^{-1}+p\times (p\ \text{mod}\ i)^{-1})(\text{mod}\ p)\\ 为什么可以加p\times (p\ \text{mod}\ i)^{-1}\\ 因为在模p的意义下p\times (p\ \text{mod}\ i)^{-1}为0\\ \therefore i^{-1}=(p-\lfloor \frac{p}{i} \rfloor) \times (p\ \text{mod}\ i)^{-1}(\text{mod}\ p) $$

用代码写成递推表达式：inv[i] = (p - p / i) * inv[p % i] % p; 初始值：inv[0] = inv[1] = 1;

LL inv[N];

void mod_inverse(LL n, LL p)
{
    inv[0] = inv[1] = 1;
    for (int i = 2; i <= n; i++)
        inv[i] = (p - p / i) * inv[p % i] % p;
}

4、中国剩余定理

可求解如下形式的一元线性同余方程组（其中 $c_1,c_2,...,c_k$ 两两互质）

$$\begin{cases} x \equiv r_1(\text{mod}\ c_1) \\ x \equiv r_2(\text{mod}\ c_2) \\ \ \ \ \ \vdots\quad \\ x \equiv r_k(\text{mod}\ c_k) \end{cases} $$

过程：

1. 计算所有模数的积 $n$。

2. 对于第 $i$ 个方程：

   a. 计算 $m_i = \frac{n}{c_i}$;

   b. 计算 $m_i$ 在模 $c_i$ 意义下的逆元 $m_i^{-1}$;

   c. 计算 $t_i=m_im_i^{-1}$（不对 $n_i$ 取模）

3. 方程组在模 $n$ 意义下的唯一解为：$x=\sum\limits^k_{i=1}a_it_i(\text{mod}\ n)$

i64 CRT(int k, std::vector<i64>&r, std::vector<i64>&c) {
    i64 n = 1, ans = 0;
    for (int i = 1; i <= k; i++) n = n * c[i];
    for (int i = 1; i <= k; i++) {
        i64 m = n / c[i], b, y;
        exgcd(m, c[i], b, y); // m * b = 1(mod c[i]);  => m*b + c[i]*y = 1
        ans = (ans + r[i] * m * b % n) % n;
    }
    return (ans % n + n) % n;
}

5、扩展中国剩余定理

可求解如下形式的一元线性同余方程组（其中 $n_1,n_2,...,n_k$ 两两不一定互质）

$$\begin{cases} x \equiv r_1(\text{mod}\ c_1) \\ x \equiv r_2(\text{mod}\ c_2) \\ \ \ \ \ \vdots\quad \\ x \equiv r_k(\text{mod}\ c_k) \end{cases} $$

前两个方程：$x\equiv r_1(\text{mod} \ c_1),x\equiv r_2(\text{mod} \ c_2)$。

转化为不定方程：$x=c_1p+r_1,x=c_2q+r_2$。

则 $c_1p-c_2q=r_2-r_1$

由裴蜀定理：

当 $\gcd(c_1,c_2) \nmid(r_2-r_1)$ 时，无解

当 $\gcd(c_1,c_2)\mid(r_2-r_1)$ 时，有解

由 $\text{exgcd}$ 算法：

特解为 $p=p\times\cfrac{r_2-r_1}{\gcd},q=q\times\cfrac{r_2-r_1}{\gcd}$

通解为 $P=p+\cfrac{c_2}{\gcd}k,Q=q-\cfrac{c_1}{\gcd}k$

所以 $x=c_1P+r_1=\cfrac{c_1c_2}{\gcd}k+c_1p+r_1$

前两个方程等价合并为一个方程 $x\equiv r(\text{mod} \ c)$

其中 $r=c_1p+r_1,c=\text{lcm}(c_1,c_2)$

所以合并 $n-1$ 次同余方程就可以求解了。

i64 EXCRT(int k, std::vector<i64>&r, std::vector<i64>&c) {
    i64 c1, c2, r1, r2, p, q;
    c1 = c[1], r1 = r[1];
    for (int i = 2; i <= k; i++) {
        c2 = c[i], r2 = r[i];
        i64 d = exgcd(c1, c2, p, q);
        if ((r2 -r1) % d) {
            return -1;
        }
        p = p * (r2 - r1) / d;
        p = (p % (c2 / d) + (c2 / d)) % (c2 / d);
        r1 = c1 * p + r1;
        c1 = c1 * c2 / d;
    }
    return (r1 * c1 + c1) % c1;
}

6、BSGS算法

给定整数 $a,b,p$，$a,p$ 互质。

求满足 $a^x \equiv b(\text{mod}\ p)$ 的最小非负整数 $x$。

时间复杂度 $O(\sqrt p)$

由扩展欧拉定理 $a^x \equiv a^{x\ \text{mod} \ \varphi(p)}(\text{mod} \ p)$

可知 $a^x$ 模 $p$ 意义下的最小循环节为 $\varphi(p)$

因为 $\varphi(p) \lt p$，所以我们考虑 $x\in [0,p]$，这样就必然能找到最小满足条件的整数 $x$

令 $x=im-j$，其中 $m=\lceil \sqrt p \rceil,i \in [1,m],j\in [0,m-1]$

则 $a^{im-j} \equiv b(\text{mod}\ p)$

即 $(a^m)^i\equiv ba^j(\text{mod}\ p)$

1. 先枚举 $j$，把 $(ba^j,j)$ 丢进哈希表，如果 $key$ 重复，就用更大的 $j$ 替换旧的
2. 再枚举 $i$，计算 $(a^m)^i$，到哈希表中查找是否有相等的 $key$，找到第一即结束。则最小的 $x=im-j$

i64 bsgs(i64 a, i64 b, i64 p) {
    a %= p, b %= p;
    if (b == 1) return 0;
    i64 m = std::ceil(std::sqrt(p));
    i64 t = b;
    std::unordered_map<i64, i64> mp;
    mp[b] = 0;
    for (int j = 1; j < m; j++) {
        t = t * a % p;
        mp[t] = j;
    }
    i64 mi = 1;
    for (int i = 1; i <= m; i++) {
        mi = mi * a % p;
    }
    t = 1;
    for (int i = 1; i <= m; i++) {
        t = t * mi % p;
        if (mp.count(t)) {
            return i * m - mp[t];
        }
    }
    return -1;
}

7、扩展BSGS算法

给定整数 $a,b,p$。

求满足 $a^x \equiv b(\text{mod}\ p)$ 的最小非负整数 $x$。

• 当 $a,p$ 互质，直接使用 BSGS 算法
• 当 $a,p$ 不互质，我们要想办法让他们变得互质

原方程可以等价为 $aa^{x-1}+py=b$

设 $d_1=\gcd(a,p)$。如果 $d_1\nmid b$，则原方程无解。（裴蜀定理）

否则方程两边同时除以 $d_1$，得到 $\cfrac{a}{d_1}a^{x-1}\equiv \cfrac{b}{d_1}\Big(\text{mod}\ \cfrac{p}{d_1}\Big)$

如果 $a$ 和 $\cfrac{p}{d_1}$ 仍不互质就再除。

设 $d_2=\gcd(a,\cfrac{p}{d_1})$。如果 $d_2 \nmid \cfrac{b}{d_1}$，则原方程无解。（裴蜀定理）

否则方程两边同时除以 $d_2$，得到 $\cfrac{a^2}{d_1d_2}a^{x-2}\equiv \cfrac{b}{d_1d_2}\Big(\text{mod}\ \cfrac{p}{d_1d_2}\Big)$

不停的判断下去直到 $a\perp \cfrac{p}{d_1 \cdots d_k}$，记 $D=\prod\limits_{i=1}^kd_i$

原方程变为 $\cfrac{a^k}{D}a^{x-k}\equiv \cfrac{b}{D}\left(\text{mod}\ \cfrac{p}{D}\right)$

因为 $a\perp \cfrac{p}{D}$，则 $\cfrac{a^k}{D}\perp\cfrac{p}{D}$，则 $\cfrac{a^k}{D}$ 就有逆元了，把它丢到方程右边，这就是一个 BSGS问题了，求解 $x-k$ 后再加上 $k$ 就是答案了。

即求解 $Aa^{im-j}\equiv b^{\prime}(\text{mod}\ p^{\prime})$，其中 $A=\cfrac{a^k}{D},b^{\prime}=\cfrac{b}{D},p^{\prime}=\cfrac{p}{D}$

i64 exbsgs(i64 a, i64 b, i64 p) {
    a %= p, b %= p;
    if (b == 1 || p == 1) return 0;
    i64 d, k = 0, A = 1;
    while (true) {
        d = std::gcd(a, p);
        if (d == 1) break;
        if (b % d) {
            return -1;
        }
        k++;
        b /= d;
        p /= d;
        A = A * (a / d) % p;
        if (A == b) {
            return k;
        }
    }
    i64 m = std::ceil(std::sqrt(p));
    i64 t = b;
    std::unordered_map<i64, i64> mp;
    mp[b] = 0;
    for (int j = 1; j < m; j++) {
        t = t * a % p;
        mp[t] = j;
    }
    i64 mi = 1;
    for (int i = 1; i <= m; i++) {
        mi = mi * a % p;
    }
    t = A;
    for (int i = 1; i <= m; i++) {
        t = t * mi % p;
        if (mp.count(t)) {
            return i * m - mp[t] + k;
        }
    }
    return -1;
}

三、数论函数

1、莫比乌斯反演

莫比乌斯函数：$\mu(x)$

定义：$\mu(x)=\begin{cases} 1 & x\ 含有平方因子 \\ (-1)^k &k 为\ x\ 的本质不同质因子个数 \end{cases}$

$n$ 的所有约数的莫比乌斯函数的和：$S(n)=\sum\limits_{d|n}\mu(d)$

性质：$S(n)=[n=1]$ 即：$\sum\limits_{d|n}\mu(d)=[n=1]$

应用：$\sum\limits_{d\mid \gcd(i,j)}\mu(d)=[\gcd(i,j)=1]$

证明：

根据整数的唯一分解定理：$n=p_1^{\alpha_1}p_2^{\alpha_2}...p_k^{\alpha_k}$，因为 $n$ 不等于 $1$，所以 $k\ge1$

$n$ 的每一个约数设为 $d$，分解 $d$ 可以得到 $d=p_1^{\beta_1}p_2^{\beta_2}...p_k^{\beta_k}，0\le \beta_i \le \alpha_i$。

我们只需要求证 $\sum \mu(d)=0$ 就是了。

由于莫比乌斯函数的性质：$d$ 若含有平方的因子，直接是零，不用加。我们就考虑另一种情况。

$\beta_i$ 取 $0$ 还是 $1$ ，我们就想到了组合数。

我们可以在 $k$ 个里取 $0$ 个 $1$ ：$C_k^0(-1)^0$

我们可以在 $k$ 个里取 $1$ 个 $1$ ：$C_k^1(-1)^1$

我们可以在 $k$ 个里取 $2$ 个 $1$ ：$C_k^2(-1)^2$

一直到我们在 $k$ 个里取 $k$ 个 $1$：$C_k^k(-1)^k$

这样我们加起来就是：$\sum\limits^k_{i=0}C_k^i(-1)^i$，由二项式定理：$0=(1-1)^k=\sum\limits^k_{i=0}C_k^i(1)^0(-1)^i=\sum\limits^k_{i=0}C_k^i(-1)^i$

证毕。

莫比乌斯函数和欧拉函数的关系

欧拉函数：对于正整数 $n$ ，欧拉函数是小于等于 $n$ 的正整数中与 $n$ 互质的数的数目，记作 $\varphi(n)$ 。特别的 $\varphi(1)=1$

即：$\varphi(n)=\sum\limits^n_{i=1}[\gcd(i,n)=1]$

关系：$\sum\limits_{d|n}\mu(d)\cfrac nd=\varphi(n)$

证明：

当 $n=1$ 时，$d=1$，$\varphi(n)=\mu(n)=1$

当 $n\gt 1$ 时，$n=p_1^{\alpha_1}p_2^{\alpha_2}...p_k^{\alpha_k}$，我们只考虑 $\alpha_i\le1$ 的情况，即 $n'=p_1p_2...p_k$

$\sum\limits_{d|n}\mu(d)\cfrac nd=n\sum\limits_{d|n'}\cfrac {\mu(d)}{d}$

知道的，约数是由质因子或乘积构成的，由容斥原理

$n\sum\limits_{d|n'}\cfrac {\mu(d)}{d}=n\left(1-\left(\cfrac 1{p_1}+...+\cfrac 1{p_k}\right)+\left(\cfrac 1{p_1p_2}+...+\cfrac 1{p_{k-1}p_k}\right)-...\right)=n\left(1-\cfrac{1}{p_1}\right)\left(1-\cfrac{1}{p_2}\right)...\left(1-\cfrac{1}{p_k}\right)=\varphi(n)$

证毕。

莫比乌斯反演

前置知识：迪利克雷卷积 *

定义：$f(n),g(n)$ 是两个积性函数，$(f*g)(n)=\sum\limits_{d|n}f(d)g\left(\cfrac nd\right)=\sum\limits_{d|n}f\left(\cfrac nd\right)g(d)$

规律：

1、交换律：$f*g=g*f$

2、结合律：$(f*g)*h=f*(g*h)$

3、分配律：$(f+g)*h=f*h+g*h$

三个常用函数：

1、元函数 $\varepsilon(n)=[n=1]$

2、常数函数 $1(n)=1$

3、恒等函数 $id(n)=n$

常用卷积关系

1、$\sum\limits_{n|d}\mu(d)=[n=1]\Rightarrow\mu*1=\varepsilon$

2、$\sum\limits_{d|n}\mu(d)=n\Rightarrow\varphi*1=id$

3、$\sum\limits_{d|n}\mu(d)\cfrac nd=\varphi(n)\Rightarrow\mu*id=\varphi$

4、$f*\varepsilon=f$

5、$f*1\neq f$

定义：若 $f(n)=\sum\limits_{d|n}g(d)$ 则 $g(n)=\sum\limits_{d|n}\mu(d)f(\cfrac n d)$

$f(n),g(n)$ 均为积性函数。

$f(n)$ 称为 $g(n)$ 的莫比乌斯变换

$g(n)$ 称为 $f(n)$ 的莫比乌斯逆变换

证明：

$$\begin{split} \sum_{d|n}\mu(d)f(\cfrac nd)&=\sum_{n|d}\mu(d)\sum_{k\mid\tfrac nd}g(k) \\ &=\sum_{d|n}\sum_{k|\tfrac nd}\mu(d)g(k)=\sum_{k|n}\sum_{d|\tfrac nk}\mu(d)g(k) \\ &=\sum_{k|n}g(k)\sum_{d|\tfrac nk}\mu(d)\\ &=g(n) \end{split} $$

可以写成另一种形式：若 $f(n)=\sum\limits_{n|d}g(d)$ 则 $g(n)=\sum\limits_{n|d}\mu(\cfrac dn)f(d)$

常见题型

P3455 ZAP-Queries：$\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k]$

$$\begin{split} &\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k] \\ =&\sum\limits_{i=1}^{\lfloor \tfrac n k\rfloor}\sum\limits_{j=1}^{\lfloor \tfrac m k\rfloor}[\gcd(i,j)=1] \end{split} $$

根据 $\sum\limits_{d\mid \gcd(i,j)}\mu(d)=[\gcd(i,j)=1]$ 可以得到

$$\begin{split} &=\sum\limits_{i=1}^{\lfloor \tfrac n k\rfloor}\sum\limits_{j=1}^{\lfloor \tfrac m k\rfloor}\sum\limits_{d\mid \gcd(i,j)}\mu(d) \\ &=\sum\limits_{i=1}^{\lfloor \tfrac n k\rfloor}\sum\limits_{j=1}^{\lfloor \tfrac m k\rfloor}\sum\limits_{d=1}^{\min(\lfloor \tfrac n k\rfloor,\lfloor \tfrac m k\rfloor)}[d|i][d|j]\mu(d) \end{split} $$

交换求和次序，先枚举 $d$ 可得

$$\begin{split} &=\sum\limits_{d=1}^{\min(\lfloor \tfrac n k\rfloor,\lfloor \tfrac m k\rfloor)}\mu(d)\sum\limits_{i=1}^{\lfloor \tfrac n k\rfloor}[d\ |\ i]\sum\limits_{j=1}^{\lfloor \tfrac m k\rfloor}[d\ |\ j]\\ &=\sum\limits_{d=1}^{\min(\lfloor \tfrac n k\rfloor,\lfloor \tfrac m k\rfloor)}\mu(d)\left\lfloor \cfrac n {kd}\right\rfloor\left\lfloor \cfrac m {kd}\right\rfloor \end{split} $$

显然可以使用数论分块解决这个问题。

数论分块（整数分块）：形如 $\sum\limits_{i=1}^nf(i)\left\lfloor\cfrac ni\right\rfloor$ 就能使用数论分块解决。

P2522 Problem b：$\sum\limits_{i=x}^n\sum\limits_{j=y}^m[\gcd(i,j)=k]$

分析可以知道，通过差分的性质，我们可以求这四个出来进行加减操作可以得到原式：

$$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k] \tag{1}$$ $$\sum_{i=1}^{x-1}\sum_{j=1}^m[\gcd(i,j)=k] \tag{2}$$ $$\sum_{i=1}^n\sum_{j=1}^{y-1}[\gcd(i,j)=k] \tag{3}$$ $$\sum_{i=1}^{x-1}\sum_{j=1}^{y-1}[\gcd(i,j)=k] \tag{4} $$

所以我们的答案就是 $(1)-(2)-(3)+(4)$

这 $4$ 个式子的推导和上面的类似。

#2185 约数个数和：$\sum\limits_{i=1}^n\sum\limits_{j=1}^md(ij)$ （$d(x)$ 为 $x$ 的约数个数）

前置知识：$d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1]$

将这个式子化简：

$$\begin{split} d(ij) &= \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1] \\ &= \sum_{x|i}\sum_{y|j}\sum_{k|\gcd(x,y)}\mu(k) \\ &= \sum_{k=1}^{\min(i,j)}\sum_{x|i}\sum_{y|j}[k|\gcd(x,y)]\mu(k) \\ &= \sum_{k=1}^{\min(i,j)}\mu(k)\sum_{x|i}\sum_{y|j}[k|\gcd(x,y)] \\ &= \sum_{k=1}^{\min(i,j)}\mu(k)\sum_{kx|i}\sum_{ky|j}[k|\gcd(kx,ky)] \\ &= \sum_{k=1}^{\min(i,j)}\mu(k)\sum_{x|\tfrac ik}\sum_{y|\tfrac jk}1 \\ &= \sum_{k=1}^{\min(i,j)}\mu(k)d\left(\frac ik\right)d\left(\frac jk\right) \end{split} $$

带入原式：

$$\begin{split} \sum_{i=1}^n\sum_{j=1}^md(ij)&=\sum_{i=1}^n\sum\limits_{j=1}^m\sum_{k=1}^{\min(i,j)}\mu(k)d\left(\frac ik\right)d\left(\frac jk\right) \\ &=\sum_{k=1}^{\min(n,m)}\sum_{i=1}^n\sum\limits_{j=1}^m[k|i][k|j]\mu(k)d\left(\frac ik\right)d\left(\frac jk\right) \\ &=\sum_{k=1}^{\min(n,m)}\sum_{i=1}^{\left\lfloor \tfrac nk \right\rfloor}\sum\limits_{j=1}^{\left\lfloor \tfrac mk \right\rfloor}\mu(k)d\left(i\right)d\left(j\right) \\ &=\sum_{k=1}^{\min(n,m)}\mu(k)\sum_{i=1}^{\left\lfloor \tfrac nk \right\rfloor}d\left(i\right)\sum\limits_{j=1}^{\left\lfloor \tfrac mk \right\rfloor}d\left(j\right) \\ \end{split} $$

我们记 $\sum\limits_{i=1}^{p}d\left(i\right)=S(p)$ 就是 $d$ 的前缀和，那么原式可以写成：

$$\sum_{i=1}^n\sum_{j=1}^md(ij)=\sum_{k=1}^{\min(n,m)}\mu(k)S\left({\left\lfloor \frac nk \right\rfloor}\right)S\left({\left\lfloor \frac mk \right\rfloor}\right) \\ $$

所以我们只需要预处理 $\mu,d$ 的前缀和，再处理分块就行了。

SPOJ LCMSUM：$\sum\limits_{i=1}^n \text{lcm}(i,n)$

化简原式得：

$$\sum_{i=1}^n \frac {in} {\gcd(i,n)}$$

这个式子可以这样化简，将原式复制一份并且颠倒顺序，将 $n$ 单独提出来（如果不懂可以自己用简单的式子推推）

$$\frac12\left(\sum_{i=1}^{n-1} \frac {in} {\gcd(i,n)}+\sum_{i=n-1}^{1} \frac {in} {\gcd(i,n)}\right)+n $$

根据 $\gcd(i,n)=\gcd(n-i,n)$，原式变为

$$\frac12\left(\sum_{i=1}^{n-1} \frac {in} {\gcd(i,n)}+\sum_{i=n-1}^{1} \frac {in} {\gcd(n-i,n)}\right)+n $$

现在看两个求和式子中，分母已经同项了，我们可以合并（自己可以带 $i=1$ 去验算）

$$\frac12\sum_{i=1}^{n-1} \frac {n^2} {\gcd(i,n)}+n$$

由于（$\gcd(n,n)=n$）我们把 $n$ 重新算上，然后就化成这样了。

$$\frac12\sum_{i=1}^{n} \frac {n^2} {\gcd(i,n)}+\frac n2 $$

因为我们知道 $\gcd(i,n)$ 有些是相同的，所以只需要统计 $\gcd(i,n)=d$ 的个数。又当 $\gcd(i,n)=d$ 时，$\gcd(\cfrac id,\cfrac nd)=1$，所以 $\gcd(i,n)=d$ 的个数一共有 $\varphi(\cfrac nd)$ 个。相当于 $\sum\limits_{\tfrac id=1}^{\tfrac nd}[\gcd(\tfrac id,\tfrac nd)=1]=\varphi(\cfrac nd)$ 欧拉函数的公式，我们每次枚举 $d$ 看每个 $d$ 加了多少次

$$\frac12\sum_{d|n} \frac {n^2\varphi(\frac nd)}{d} +\frac n2 $$

改变求和次序，并且令 $d'=\cfrac nd$ 合并公式

$$\frac12n\left(\sum_{d'|n} {d'\varphi(d')}+1\right) $$

设 $\text{g}(n)=\sum\limits_{d|n}d\cdot \varphi(d)$，已知 $\text{g}$ 为积性函数，可以 $O(n)$ 的时间筛出。每次询问 $O(1)$ 计算。

下面是推导 $\text{g}$ 的过程，考虑 $\text{g}(p_j^k)$ 它的约数只有 $p_j^0,p_j^1,...,p_j^k$，因此

$$\text{g}(p_j^k)=\sum_{w=1}^kp_j^w \cdot \varphi(p_j^w) $$

又因为 $\varphi(p_j^{w})=p_j^{w} - p_j^{w-1}=p_j^{w}\cdot (1- \frac{1}{p_j})=p_j^{w-1}\cdot (p_j-1)$，原式可以化简

$$\sum_{w=0}^kp_j^{2w-1}\cdot(p_j-1)$$

于是有

$$\text{g}(p_j^{k+1})=\text{g}(p_j^{k})+p_j^{2k+1}\cdot(p_j-1) $$

对于 $i$ 是质数的时候

$$\begin{split} \text{g}(i)=&1\cdot \varphi(1)+i\cdot \varphi(i) \\ =&i\cdot (i-1) + 1 \end{split} $$

现在对于 $p_j|i$ 的情况，这个 $i$ 就是欧拉筛的那个 $i$，令 $i=a\cdot p_j^w(\gcd(a,p_j)=1)$

$$\text{g}(i\cdot p_j)=\text{g}(a)\cdot \text{g}(p_j^{w+1}) \\ \text{g}(i)=\text{g}(a) \cdot \text{g}(p_j^w) $$

即

$$\text{g}(i\cdot p_j)-\text{g}(i)=\text{g}(a)\cdot p_j^{2w+1} \cdot (p_j-1) $$

同理有

$$\text{g}(i)-\text{g}(\frac i {p_j})=\text{g}(a)\cdot p_j^{2w-1} \cdot (p_j-1) $$

因此

$$\text{g}(i\cdot p_j)=\text{g}(i)+\left(\text{g}(i)-\text{g}(\frac i {p_j})\right)\cdot p_j^2 $$

P1829 Crash的数字表格：$\sum\limits_{i=1}^n\sum\limits_{j=1}^m\text{lcm}(i,j)(\text{mod}\ 20101009)$

和上题有出入，首先先化简

$$\begin{split} \sum_{i=1}^n\sum_{j=1}^m\text{lcm}(i,j)=&\sum_{i=1}^n\sum_{j=1}^m\frac {i\cdot j}{\gcd(i,j)}\\ =&\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^n\frac{i\cdot j}{d}[\gcd(i,j)=d] \\ =&\sum_{id=1}^n\sum_{jd=1}^m\sum_{d=1}^n\frac{id\cdot jd}{d}[\gcd(id,jd)=d] \\ =&\sum_{i=1}^{\left\lfloor\tfrac nd\right\rfloor}\sum_{j=1}^{\left\lfloor\tfrac md\right\rfloor}\sum_{d=1}^nd\cdot i\cdot j[\gcd(i,j)=1] \\ =&\sum_{d=1}^nd\sum_{i=1}^{\left\lfloor\tfrac nd\right\rfloor}\sum_{j=1}^{\left\lfloor\tfrac md\right\rfloor}i\cdot j[\gcd(i,j)=1]\\ \end{split} $$

我们记 $S(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}i\cdot j[\gcd(i,j)=1]$ 化简它

$$\begin{split} S(n,m)=&\sum_{i=1}^n\sum_{j=1}^m\sum_{d|\gcd(i,j)}\mu(d)\cdot i\cdot j \\ =&\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^n[d|i][d|j]\mu(d)\cdot i\cdot j \\ =&\sum_{d=1}^n\mu(d)\sum_{id=1}^n\sum_{jd=1}^mid\cdot jd \\ =&\sum_{d=1}^n\mu(d)\cdot d^2\sum_{i=1}^{\left\lfloor\tfrac nd\right\rfloor}\sum_{j=1}^{\left\lfloor\tfrac md\right\rfloor}i\cdot j \end{split} $$

所以我们观察，前面一部分可以通过数论分块求出来，后面的可以直接算出来，定义一个函数 $f(n,m)=\sum\limits_{i=1}^n\sum\limits_{j=1}^mi\cdot j$

$$f(n,m)=\sum\limits_{i=1}^n\sum\limits_{j=1}^mi\cdot j=\frac{n\cdot (n+1)}{2} \times \frac{m\cdot (m+1)}{2} $$

所以 $S(n,m)$ 函数化简的结果为

$$S(n,m)=\sum_{d=1}^n\mu(d)\cdot d^2\cdot f\left({\left\lfloor\frac nd\right\rfloor},{\left\lfloor\frac md\right\rfloor}\right) $$

将 $S(n,m)$ 带回原式得：

$$\sum_{d=1}^nd\cdot S\left({\left\lfloor\frac nd\right\rfloor},{\left\lfloor\frac md\right\rfloor}\right) $$

这个式子我们又可以通过数论分块得到。上面的推导都基于 $n\le m$。

时间复杂度为 $O(n+m)$ 瓶颈为线性筛