算术和算法

正整数

从Al-Khwarizmi开始，算术（Arithmetic）分为两个分支
- 代数（Algebra），之后发展为分析（Analysis），它对应于数学
- 算法（Algorithm），之后发展为AI（Artificial Intelligence），它对应于计算机
在这里，我们考虑算术、算法。最初，算法源于算术，算术从正整数开始
$\mathbb{Z}_+ = \{ 1, 2, 3, \ldots \}.$
如果对于正整数、、，我们有
$n = ab,$
那么我们称、整除（Divide），记为、
- $a$ 、 $b$ 称为 $n$ 的除子（Divisor）或因子（Factor）。在代数几何中，我们通常使用除子；在算术中，我们通常使用因子
- $n$ 称为 $a$ 、 $b$ 的倍数（Multiple）
- 如果正整数 $p$ 的因子只有 $1$ 、 $p$ ，那么 $p$ 称为素数（Prime Number）

import math

// (n) -> (factor)
def find_factor(n):
    factor = 1

    // If n = ab, 1 < a <= b < n, then 2 <= a <= floor(sqrt(n))
    upper_bound = math.floor(math.sqrt(n))

    for i in range(2, upper_bound + 1):
        if n % i == 0:
            factor = i
            break
    
    return factor

// (n) -> (result)
def is_prime(n):
    result = False

    if find_factor(n) == 1:
        result = True

    return result

print(is_prime(23))
print(is_prime(233))
print(is_prime(2333))
print(is_prime(23333))
print(is_prime(233333))
print(find_factor(233333))

// ----- Run Program -----

> python arithmetic.py
True
True
True
True
False
353

Euclid算法

现在，我们考虑整数，即增加零、负数
$\mathbb{Z} = \{ 0, \pm 1, \pm 2, \pm 3, \ldots \}.$
如果对于整数、、，我们有
$n = ab,$
那么我们称、整除，记为、
- $a$ 、 $b$ 称为 $n$ 的除子或因子
- $n$ 称为 $a$ 、 $b$ 的倍数
- 如果整数 $p$ 的因子只有 $\pm 1$ 、 $\pm p$ ，那么 $p$ 称为素数。注意， $\pm 1$ 不是素数，而是单位（Unit），所以上面的is_prime()不适用于 $n = 1$
（带余除法）设、为整数，
- （商）考虑如下序列
  $a, a \pm b, a \pm 2b, a \pm 3b, \ldots$
  存在一个最小的非负整数
  $r = a + qb,\; q \in \mathbb{Z},$
  其中， $q$ 称为商（Quotient）
- （余数）注意， $r - b$ 、 $r + b$ 的符号相反（否则， $r$ 不是最小的非负整数）。因此，
  $(r - b)(r + b) < 0,\; 0 \leq r < |b|,$
  其中， $r$ 称为余数（Remainder）
- 我们也称 $a$ 为被除数（Dividend）， $b$ 为除数（Divisor）。如果除数为0，那么一种方法是定义商为任意整数，余数为被除数
（Euclid算法）设、为整数
- （最大公因子）如果正整数 $d$ 满足
  $d | a,\; d | b,$
  那么， $d$ 称为 $a$ 、 $b$ 的公因子，并且最大公因子（Greatest Common Divisor，GCD）记为 $gcd(a, b)$
- （最小公倍数）如果正整数 $m$ 满足
  $a | m,\; b | m,$
  那么， $m$ 称为 $a$ 、 $b$ 的公倍数，并且最小公倍数（Least Common Multiple，LCM）记为 $lcm(a, b)$
- 现在，我们希望寻找 $gcd(a, b)$ 。利用带余除法，
  $gcd(a, b) = gcd(a - qb, b) = gcd(b, r),$
  其中， $0 \leq r < |b|$ 。再次利用带余除法，
  $gcd(b, r) = gcd(b - q'r, r) = gcd(r, r'),$
  其中， $0 \leq r' < r$ 。因此，除数会递减至0，并且剩下的被除数会等于 $gcd(a, b)$
（Bézout定理）进一步，在Euclid算法中，所有项（包括）都具有的形式。因此，我们可以得到Bézout定理
$gcd(a, b) = ma + nb \text{ for some } m, n \in \mathbb{Z}.$
并且，我们可以同时计算、、
- 令
  $a_1 = a,\; b_1 = b.$
  利用带余除法，对于 $i \geq 1$ ，令
  $a_{i + 1} = b_i,\; b_{i + 1} = r_i, \text{ where } a_i = q_ib_i + r_i.$
- 被除数、除数
  $\begin{equation*}\begin{split} a_i &= m_i \cdot a + n_i \cdot b, \\ b_i &= o_i \cdot a + p_i \cdot b. \end{split}\end{equation*}$
- 下一个被除数、下一个除数
  $\begin{equation*}\begin{split} a_{i + 1} = b_i &\Rightarrow m_{i + 1} = o_i,\; n_{i + 1} = p_i, \\ b_{i + 1} = r_i &\Rightarrow o_{i + 1} = m_i - q_io_i,\; n_{i + 1} = n_i - q_ip_i. \end{split}\end{equation*}$

// (a, b) -> (q, r)
def division_with_remainder(a, b):
    // Python does not behave as the above definition, so we only consider positive integers
    q = a // b
    r = a % b

    return (q, r)

// (a, b) -> (d, m, n)
def Euclidean_algorithm(a, b):
    d = 0; m = 0; n = 0
    
    dividend = a; divisor = b

    // dividend = mi * a + ni * b
    mi = 1; ni = 0

    // divisor = oi * a + pi * b
    oi = 0; pi = 1

    while(divisor != 0):
        (qi, ri) = division_with_remainder(dividend, divisor)
        
        dividend = divisor
        divisor = ri

        (mi, ni, oi, pi) = (oi, pi, mi - qi * oi, ni - qi * pi)

    d = dividend
    m = mi
    n = ni

    return (d, m, n)

print(division_with_remainder(2333, 233))
print(division_with_remainder(233, 2333))

a = 233; b = 2333
(d, m, n) = Euclidean_algorithm(a, b)
print(d, m, n)
print(m * a + n * b)

// ----- Run Program -----

> python arithmetic.py
(10, 3)
(0, 233)
1 -781 78
1

素因子分解

（存在性）对任意整数，我们有素因子分解（Prime Factorization）
$n = u \cdot p_1 \cdots p_k,$
其中，为单位，为素数
- 如果 $n$ 是正整数（或者负整数），那么 $u = +1$ （或者 $u = -1$ ）。因此，从下面开始，我们通常只考虑正整数的情形
- 如果 $n$ 是素数，那么它不能继续进行分解。如果 $n$ 不是素数，那么
  $n = ab,\; 1 < a \leq b < n,$
  并且 $a$ 、 $b$ 可以递归地进行分解，直到不能继续进行分解
素因子分解的存在性有如下推论
- 素数有无限多个
  - 若不然，则令 $p_1, \ldots, p_k$ 为所有素数。考虑
    $n = p_1 \cdots p_k + 1.$
  - 利用素因子分解，
    $n = p_1 \cdots p_k + 1 = p_{i_1} \cdots p_{i_l}.$
    因此， $p_{i_1} | 1, \cdots, p_{i_l} | 1$ ，矛盾
- 我们可以收集重复的素因子，得到幂的形式
  $n = u \cdot \prod_i p_i^{k_i}.$
  - 如果正整数 $a$ 、 $b$ 用相同的素因子来表示（缺少的素因子指数为0），
    $a = \prod_i p_i^{k_i},\; b = \prod_i p_i^{l_i},$
    那么，对于 $m_i = \min\{ k_i, l_i \}$ 、 $n_i = \max\{ k_i, l_i \}$ ，我们有
    $gcd(a, b) = \prod_i p_i^{m_i},\; lcm(a, b) = \prod_i p_i^{n_i}.$
  - 在利用Euclid算法找到 $gcd(a, b)$ 之后，我们可以利用如下关系找到 $lcm(a, b)$
    $m_i + n_i = k_i + l_i \Rightarrow gcd(a, b) \cdot lcm(a, b) = ab.$

// (n) -> (prime_factors)
def prime_factorization(n):
    prime_factors = []

    recursive_factor(n, prime_factors)

    return prime_factors

// (n, prime_factors) -> ()
def recursive_factor(n, prime_factors):
    if is_prime(n):
        prime_factors.append(n)
    else:
        a = find_factor(n)
        b = n // a
        recursive_factor(a, prime_factors)
        recursive_factor(b, prime_factors)

print(prime_factorization(233333))
print(prime_factorization(123456789))

// ----- Run Program -----

> python arithmetic.py
[353, 661]
[3, 3, 3607, 3803]

Eratosthenes筛法

尽管素数有无限多个，然而我们可以寻找不大于给定正整数的所有素数。一种方法是Eratosthenes筛法（Sieve）
- 创建一个从 $2$ 到 $n$ 的列表
- 令 $p = 2$ 。从列表中去掉 $2p, 3p, \ldots, \lfloor{n / p}\rfloor \cdot p$
- 在列表中，找到大于 $p$ 的最小正整数。将它赋值给 $p$ ，并且重复上述过程
- 在列表中，如果不存在大于 $p$ 的正整数，那么列表包含了不大于 $n$ 的所有素数
在解析数论中，筛法可以应用于素数分布的问题（比如Goldbach猜想的较弱版本）。关于素数分布，一个经典结果是素数定理
$\pi(n) \sim \frac{n}{\ln n} \text{ as } n \to \infty,$
其中， $\pi(n)$ 是不大于 $n$ 的所有素数的个数

// (n) -> (all_primes)
def Eratosthenes_sieve(n):
    // Create a list from 2 to n
    all_primes = list(range(2, n + 1))

    p = 2

    while(max(all_primes) > p):
        // Remove 2p, 3p, ..., floor(n / p) * p from the list
        for i in range(2, math.floor(n / p) + 1):
            if i * p in all_primes:
                all_primes.remove(i * p)
        
        // Find the smallest number greater than p in the list
        for i in range(len(all_primes)):
            if all_primes[i] > p:
                p = all_primes[i]
                break

    return all_primes

print(Eratosthenes_sieve(233))

n = 233
print(len(Eratosthenes_sieve(n)) / (n / math.log(n)))
n = 2333
print(len(Eratosthenes_sieve(n)) / (n / math.log(n)))
n = 23333
print(len(Eratosthenes_sieve(n)) / (n / math.log(n)))

// ----- Run Program -----

> python arithmetic.py
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233]
1.1931457559306897
1.146782702034888
1.1215847730217052