ACM演算法——數學專題

高明y發表於2024-08-19

一、素數

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\) 的字首和,再處理分塊就行了。

相關文章