Mathematics in ACM

素数與歐拉函數

素性判定 & 素因子分解

给出一个数 n，如何判断这个数是否是素数？

$O(\sqrt{n})$

boolean isPrime(int n)
{
    if (n < 2)
        return false;
    for (int i = 2; i <= n / i; i++)
        if (n % i == 0)
            return false;
    return true;
}

将一个数 n 素因子分解，$O(\sqrt{n})$：

ArrayList<Integer> factor = new ArrayList<Integer>();
ArrayList<Integer> power = new ArrayList<Integer>();

void makeFactors(int n, ArrayList<Integer> factor, ArrayList<Integer> power)
{
    int factorCnt = 0;
    for (int i = 2; i <= n / i; ++i) {
        if ((n % i) != 0)
            continue;
        factor.add(i);
        power.add(0);
        for (; n % i == 0; n /= i)
            power.set(factorCnt, power.get(factorCnt) + 1);
        factorCnt++;
    }
    if (n > 1) {
        factor.add(n);
        power.add(1);
        factorCnt++;
    }
    // for (int i = 0; i < factorCnt; ++i)
    //    System.out.println(factor.get(i) + " " + power.get(i));
}

問題：下面代碼的複雜度是多少？

for (int i = 1; i <= n; ++i)
    makeFactors(i, factors[i], powers[i]);

篩法

埃拉托斯特尼篩法。

final int N = 123456;
final int P = 12345;
boolean[] isPrime = new boolean[N];
int[] prime = new int[P];

void makePrime()
{
    isPrime[0] = isPrime[1] = false;
    for (int i = 2; i < N; ++i)
        isPrime[i] = true;
    prime[0] = 0; // 用此數記錄質數數目
    for (int i = 2; i < N; ++i) {
        if (isPrime[i])
            prime[++prime[0]] = i;
        for (int j = i << 1; j < N; j += i)
            isPrime[j] = false;
    }
}

  0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
      2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
2)    2 3   5   7   9    11    13    15    17    19    21    23    25
3)    2 3   5   7        11    13          17    19          23    25
4)
5)    2 3   5   7        11    13          17    19          23
6)

算法的時間複雜度是多少？

$$\lim_{n\to \infty}\sum\limits_{k = 1}^{n}\frac{1}{k} - \ln(n) = \gamma\\ \gamma \approx 0.57721566490153286060651209$$

優化 $O(n)$：

void makePrime()
{
    isPrime[0] = isPrime[1] = false;
    for (int i = 2; i < N; ++i)
        isPrime[i] = true;
    prime[0] = 0; // 用此數記錄質數數目
    for (int i = 2; i < N; ++i) {
        if (isPrime[i])
            prime[++prime[0]] = i;
        for (int j = 1; j <= prime[0] && i * prime[j] < N; ++j) {
            isPrime[i * prime[j]] = false;
            if (i % prime[j] == 0)
                break;
        }
    }
}

得到素數表之後，對一個數進行素因數分解，就變得很快。

歐拉函數

定義： $\Phi(n) = |\{x: 1 \le x \le n, \gcd(n,x) = 1\}|$。

如果 $p$ 是素數，那麼 $\Phi(p) = p-1$，$\Phi(p^k) = p^k - p^{k-1}$。

$p_1,p_2,...$ 設為素數。

$\Phi(p_1^{k_1}p_2^{k_2})=(p_1^{k_1} - p_1^{k_1-1})(p_2^{k_2} - p_2^{k_2-1})$

$\Phi(\Pi p_i^{k_i}) = \Pi(p_i^{k_i} - p_i^{k_i - 1}) = n\Pi(1-\frac{1}{p_i})$

求歐拉函數/歐拉函數表。

1. 篩選素數之後素因數分解
2. 篩法

final int N = 12345;
int[] euler = new int[N];

void eulerTable()
{
    for (int i = 0; i < N; ++i)
        euler[i] = i;
    for (int i = 2; i < N; ++i) {
        if (euler[i] != i)
            continue;
        for (int j = i; j < N; j += i)
            euler[j] -= euler[j] / i;
    }
}

歐拉函數的應用在後面模運算部分會有介紹。

模運算

模集合

定義：對於加減運算自我封閉的數集。形式話的說，如果 $M$ 是模，那麼，

$$\forall a, b \in M \implies a \pm b \in M$$

例

1. 空集是模，僅有一個 0 的集合是模（但因為太平凡我們不考慮這些狀況）。
2. 整數集合是模，因為任意兩個數的差或者和都是整數。
3. 偶數集合是模，因為任意兩個偶數的差或者和都是偶數。
4. 奇數集合不是模，因為 3 - 1 = 2 不是奇數。
5. $x \in Z, M(x) = \{kx:k \in Z\}$ 是模。

這裡我們只考慮整數模。

模的性質

1. 任何模中有 0。
2. 設 $M$ 是模，$a, b \in M \implies an+bm \in M, n, m \in Z$。
3. $G(a, b) = \{an + bm: n, m \in Z\}$ 是模。
4. 任何模都是其中最小正整數的全部倍數集合。（考查 $r = n - d \lfloor \frac{n}{d} \rfloor$）。

定義：對於兩個整數 $a, b$，模 $\{an+bm:n, m \in Z\}$ 中的最小正整數稱作 $a, b$ 的最大公約數，用 $(a, b)$ 表示。

$(a, b)$ 的性質

1. $\exists x, y \in Z, s.t. (a, b) = ax + by$
2. $(a, b)|ax+by$
3. $e|a, e|b \implies e|(a,b)$
4. $(a, b)$ 是 $a, b$ 的公共因子中最大的一個。
5. $[a, b](a, b) = ab$

歐幾里得算法

問題：求 $(a, b)$。

考查 $G(a, b) = \{an+bm:n, m \in Z\}$，回顧性質 4，如果 $b \ne 0$，

$$(a, b) = (b, a\!\!\!\!\mod b)$$

所以得到算法：

int gcd(int a, int b)
{
    if (b == 0)
        return a;
    return gcd(b, a % b);
}    

一次不定方程

對於一次不定方程 $ax+by=n$。

1. 有整數解的充要條件是 $(a,b)|n$。
2. 如果 $(a, b) = 1$，$x_0,y_0$ 是方程的解。那麼方程的全部解可以表示為： $$x=x_0+bt \quad y=y_0-at$$
3. 設 $(a, b) = 1,a > 0,b > 0$，那麼大於 $ab-a-b$ 的數可以表成 $ax+by,x\ge 0, y\ge 0$；$ab-a-b$ 不能這樣表示。

擴展歐幾里得算法

class LinearFunction
{
    int x, y;
    int g, tmp;

    int gcd(int a, int b)
    {
        if (b == 0) {
            x = 1;
            y = 0;
            return a;
        }
        g = gcd(b, a % b);
        tmp = x;
        x = y;
        y = tmp - a / b * y;
        return g;
    }
}

gcd(a, b) 使得 $x, y$ 是 $ax + by = (a, b)$ 的解。那麼 13 行之後，$x, y$ 是 $bx + (a\%b)y = (b, a\%b)$ 的解。

那麼 $bx + (a- b\lfloor \frac{a}{b} \rfloor)y=(a,b)$，$ay+(x-\lfloor \frac{a}{b} \rfloor y)b=(a,b)$

那麼 $y, x-\lfloor \frac{a}{b} \rfloor y$ 是 $ax+by=(a,b)$ 的解。

找到 $x_t,y_t$ 是 $ax+by=(a,b)$ 的解之後，它也就是 $\frac{a}{(a,b)}x+\frac{b}{(a,b)}y = 1$ 的解。利用性質 2 我們得到了後面方程的所有解。

我們關心 $ax+by=n$ 的解，當 $(a, b)|n$ 的時候，它有解，並且可以利用前面解決的方程的解通過線性變換得到所有解。

模逆元

In modular arithmetic, the modular multiplicative inverse of an integer a modulo m is an integer x such that

$$a\,x \equiv 1 \pmod{m}.$$

我們在約定 $m$ 之後，可以把 $x$ 記錄為 $a^{-1}$。

我們發現這就是要找 $ax+my=1$ 的解，我們一般使用 $1 \le x < m$ 的解 $x$ 作為模逆元。故 $(a,m) = 1$ 的時候才存在模逆元。

所以我們可以使用擴展歐幾里得算法解決模逆元。

模逆元表

現在對於一個模 $M$ （一般取素數），要處理 $1, 2, ..., N (N < M)$ 的逆元。

除了依次使用擴展歐幾里得算法以外，還可以 $O(N)$ 的計算。

發現：$M \% i = M - \lfloor \frac{M}{i} \rfloor i = -\lfloor \frac{M}{i} \rfloor i \pmod{M}$，那麼 $i^{-1}=(M\%i)^{-1}(-\lfloor \frac{M}{i} \rfloor)=(M\%i)^{-1}(M-\lfloor \frac{M}{i} \rfloor)\pmod{M}$

final int N = 12345;
long[] inv = new long[N];

void invTable(int M)
{
    inv[1] = 1;
    for (int i = 2; i < N; ++i)
        inv[i] = inv[M % i] * (M - M / i) % ((long) M);
}

求組合數的模值

$$\binom{n}{k} = \frac{n!}{k!(n-k)!}$$

設 $f_i = i!\mod M$，$g_i = (i!)^{-1}\mod M$。那麼 $\binom{n}{k} = f_n g_k g_{n-k}\mod M$。

$f_i = f_{i - 1} \times i \mod M$
$g_i = g_{i - 1} \times i^{-1} \mod M$

歐拉定理

$$(a, m) = 1 \implies a^{\Phi(m)} =1\pmod{m}$$

費馬小定理

$p$ 是素數，那麼 $a^{p-1} = 1\pmod{p}$。

指數循環節

1. $(a, b) = 1$
   $a ^ b = a^{b\mod \Phi(c)} \pmod{c}$
2. $b > \Phi(c)$
   $a ^ b = a^{b\mod \Phi(c) + \Phi(c)} \pmod{c}$

矩陣與快速冪

複習：二進制冪

long modPow(long a, long b, long M)
{
    /**@return
     * pow(a, b) % M
     */
    long result = 1L % M;
    for (; b != 0; b >>= 1) {
        if ((b & 1) != 0)
            result = result * a % MOD;
        a = a * a % MOD;
    }
    return result;
}

類似的我們可以實現二進制乘法

long modSum(long a, long b, long M)
{
    /**@return
     * (a * b) % M
     */
    long result = 0;
    for (; b != 0; b >>= 1) {
        if ((b & 1) != 0)
            result = (result + a) % M;
        a = (a << 1) % M;
    }
    return result;
}    

思考：二進制的加法有何意義？

快速冪與矩陣結合

斐波那契數列

$$\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^ n = \begin{pmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{pmatrix}$$

複雜度 $O(2^3\log(n))$

線性遞推

$f_0 = a, f_1 = b,f_i = pf_{i-1}+qf_{i-2}$，求 $f_n \% M$。

$$\begin{pmatrix} p & q \\ 1 & 0 \end{pmatrix} \begin{pmatrix} f_{n-1}\\ f_{n-2} \end{pmatrix} = \begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix}$$

$$\begin{pmatrix} p & q \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} f_1\\ f_0 \end{pmatrix} = \begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix}$$

例三

設 $s_n = f_0 + f_1 + ... + f_n$，求 $s_n \% M$。

$$\begin{pmatrix} 1 & 1 & 0\\ 1 & 0 & 0\\ 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} f_{n - 1}\\ f_{n - 2}\\ s_{n - 1}\\ \end{pmatrix} = \begin{pmatrix} f_{n}\\ f_{n - 1}\\ s_{n}\\ \end{pmatrix}$$

$$\begin{pmatrix} 1 & 1 & 0\\ 1 & 0 & 0\\ 1 & 1 & 1 \end{pmatrix}^{n-1} \begin{pmatrix} f_{1}\\ f_{0}\\ s_{1}\\ \end{pmatrix} = \begin{pmatrix} f_{n}\\ f_{n - 1}\\ s_{n}\\ \end{pmatrix}$$

作業

1. 求 $|\{p:a\le p\le b,p\ 是素數\}|$，其中 $b = 5527326628, a=b-10^5$。

2. 求 $\sum\limits_{i = 1}^{15527326628} f_i \mod (10^{12}+7)$。

作者

張靜之