跳到主要内容

数学算法

数学算法是算法竞赛和实际编程中的重要基础。本章将深入讲解数论基础、质数筛法、快速幂等核心算法，帮助你理解其数学原理和实际应用。

最大公约数与最小公倍数​

最大公约数 (GCD)​

最大公约数（Greatest Common Divisor）是指能够同时整除两个或多个整数的最大正整数。例如，gcd(12, 18) = 6，因为 6 是同时整除 12 和 18 的最大整数。

欧几里得算法原理​

欧几里得算法（又称辗转相除法）是计算 GCD 的经典方法，其核心基于以下数学定理：

gcd⁡(a,b)=gcd⁡(b,a mod b)\gcd(a, b) = \gcd(b, a \bmod b)

证明过程：

设 g=gcd⁡(a,b)g = \gcd(a, b)，则 g∣ag \mid a 且 g∣bg \mid b。

由于 a mod b=a−⌊a/b⌋⋅ba \bmod b = a - \lfloor a/b \rfloor \cdot b，而 gg 能整除 aa 和 bb，所以 gg 也能整除 a mod ba \bmod b。

反过来，如果 g∣bg \mid b 且 g∣(a mod b)g \mid (a \bmod b)，由于 a=⌊a/b⌋⋅b+(a mod b)a = \lfloor a/b \rfloor \cdot b + (a \bmod b)，所以 g∣ag \mid a。

这证明了 {a,b}\{a, b\} 和 {b,a mod b}\{b, a \bmod b\} 的公约数集合完全相同，因此它们的最大公约数也相同。

递归实现​

/**
 * 计算两个数的最大公约数（递归版本）
 * 时间复杂度: O(log(min(a, b)))
 * 
 * @param a 第一个数
 * @param b 第二个数
 * @return a 和 b 的最大公约数
 */
public int gcd(int a, int b) {
    // 基本情况：如果 b 为 0，则 a 就是 GCD
    if (b == 0) {
        return a;
    }
    // 递归调用：gcd(b, a % b)
    return gcd(b, a % b);
}

迭代实现​

/**
 * 计算两个数的最大公约数（迭代版本）
 * 避免递归栈溢出，适合处理大数
 */
public int gcdIterative(int a, int b) {
    while (b != 0) {
        int temp = a % b;
        a = b;
        b = temp;
    }
    return a;
}

时间复杂度分析​

根据拉梅定理（Lamé's Theorem），欧几里得算法的递归调用次数不超过 5⋅log⁡10(min⁡(a,b))5 \cdot \log_{10}(\min(a, b))。更精确地说：

  • 如果 a>b≥1a > b \geq 1 且 b<Fnb < F_n（FnF_n 为第 nn 个斐波那契数），则最多进行 n−2n-2 次递归调用
  • 最坏情况发生在连续的斐波那契数上，例如 gcd⁡(Fn,Fn−1)\gcd(F_n, F_{n-1})

因此，时间复杂度为 O(log⁡min⁡(a,b))O(\log \min(a, b))。

最小公倍数 (LCM)​

最小公倍数（Least Common Multiple）是能够被两个或多个整数整除的最小正整数。例如，lcm(12, 18) = 36。

数学关系​

LCM 与 GCD 存在如下关系：

lcm(a,b)=a×bgcd⁡(a,b)\text{lcm}(a, b) = \frac{a \times b}{\gcd(a, b)}

证明：

设 a=p1a1⋅p2a2⋯pkaka = p_1^{a_1} \cdot p_2^{a_2} \cdots p_k^{a_k}，b=p1b1⋅p2b2⋯pkbkb = p_1^{b_1} \cdot p_2^{b_2} \cdots p_k^{b_k}（质因数分解）

则：

  • gcd⁡(a,b)=p1min⁡(a1,b1)⋅p2min⁡(a2,b2)⋯pkmin⁡(ak,bk)\gcd(a, b) = p_1^{\min(a_1, b_1)} \cdot p_2^{\min(a_2, b_2)} \cdots p_k^{\min(a_k, b_k)}
  • lcm(a,b)=p1max⁡(a1,b1)⋅p2max⁡(a2,b2)⋯pkmax⁡(ak,bk)\text{lcm}(a, b) = p_1^{\max(a_1, b_1)} \cdot p_2^{\max(a_2, b_2)} \cdots p_k^{\max(a_k, b_k)}

由于 min⁡(ai,bi)+max⁡(ai,bi)=ai+bi\min(a_i, b_i) + \max(a_i, b_i) = a_i + b_i，所以 gcd⁡(a,b)⋅lcm(a,b)=a⋅b\gcd(a, b) \cdot \text{lcm}(a, b) = a \cdot b。

实现​

/**
 * 计算两个数的最小公倍数
 * 注意：先除后乘，避免整数溢出
 * 
 * @param a 第一个数
 * @param b 第二个数
 * @return a 和 b 的最小公倍数
 */
public long lcm(int a, int b) {
    // 先除后乘，避免 a * b 溢出
    return (long) a / gcd(a, b) * b;
}

扩展欧几里得算法​

扩展欧几里得算法不仅能求出 gcd⁡(a,b)\gcd(a, b)，还能找到整数 xx 和 yy，使得：

ax+by=gcd⁡(a,b)ax + by = \gcd(a, b)

这在求解线性同余方程和模逆元时非常有用。

/**
 * 扩展欧几里得算法
 * 返回 gcd(a, b)，同时求出 x 和 y 使得 ax + by = gcd(a, b)
 */
public int extendedGcd(int a, int b, int[] xy) {
    if (b == 0) {
        xy[0] = 1;  // x = 1
        xy[1] = 0;  // y = 0
        return a;
    }
    
    int gcd = extendedGcd(b, a % b, xy);
    // 根据递归结果更新 x 和 y
    int temp = xy[0];
    xy[0] = xy[1];
    xy[1] = temp - (a / b) * xy[1];
    
    return gcd;
}

/**
 * 求解 a 在模 m 下的乘法逆元
 * 逆元存在的条件：gcd(a, m) = 1
 */
public int modInverse(int a, int m) {
    int[] xy = new int[2];
    int gcd = extendedGcd(a, m, xy);
    if (gcd != 1) {
        return -1;  // 逆元不存在
    }
    return (xy[0] % m + m) % m;  // 确保结果为正
}

质数与筛法​

质数判断​

质数（素数）是只能被 1 和自身整除的大于 1 的正整数。

基本判断方法​

判断一个数 nn 是否为质数，只需检查 22 到 n\sqrt{n} 之间是否存在因子：

/**
 * 判断 n 是否为质数
 * 时间复杂度: O(√n)
 */
public boolean isPrime(int n) {
    // 小于 2 的数都不是质数
    if (n < 2) return false;
    
    // 2 是唯一的偶质数
    if (n == 2) return true;
    
    // 排除其他偶数
    if (n % 2 == 0) return false;
    
    // 只需检查奇数因子
    for (int i = 3; i * i <= n; i += 2) {
        if (n % i == 0) return false;
    }
    
    return true;
}

为什么只需检查到 √n？​

如果 n=a×bn = a \times b，其中 a≤ba \leq b，则必然有 a≤na \leq \sqrt{n}。因此，如果在 n\sqrt{n} 以内没有找到因子，那么 n\sqrt{n} 以外也不会有因子。

埃拉托斯特尼筛法​

埃拉托斯特尼筛法（简称埃氏筛）是一种高效找出 11 到 nn 所有质数的算法。

算法原理​

  1. 初始化一个布尔数组，标记所有数为"可能是质数"
  2. 从 2 开始，将 2 的所有倍数标记为合数
  3. 找到下一个未被标记的数，将其所有倍数标记为合数
  4. 重复直到处理完所有数

实现​

/**
 * 埃氏筛法求 1 到 n 的所有质数
 * 时间复杂度: O(n log log n)
 * 空间复杂度: O(n)
 */
public int countPrimes(int n) {
    if (n <= 2) return 0;
    
    // isPrime[i] = true 表示 i 是质数
    boolean[] isPrime = new boolean[n];
    Arrays.fill(isPrime, true);
    
    // 0 和 1 不是质数
    isPrime[0] = false;
    isPrime[1] = false;
    
    int count = 0;
    
    for (int i = 2; i < n; i++) {
        if (isPrime[i]) {
            count++;
            
            // 从 i*i 开始标记，因为 i*k (k < i) 已经被更小的质数标记过
            // 使用 long 防止 i*i 溢出
            if ((long) i * i < n) {
                for (int j = i * i; j < n; j += i) {
                    isPrime[j] = false;
                }
            }
        }
    }
    
    return count;
}

时间复杂度证明​

设 pp 为质数，内层循环执行的次数约为 n/pn/p。

总时间复杂度为：

O(∑p≤nnp)=O(n⋅∑p≤n1p)=O(nlog⁡log⁡n)O\left(\sum_{p \leq n} \frac{n}{p}\right) = O\left(n \cdot \sum_{p \leq n} \frac{1}{p}\right) = O(n \log \log n)

最后一步利用了质数倒数和的增长性质：∑p≤n1p≈log⁡log⁡n\sum_{p \leq n} \frac{1}{p} \approx \log \log n。

优化：只筛奇数​

由于除了 2 以外的偶数都不是质数，可以只处理奇数来减少空间和时间：

public int countPrimesOptimized(int n) {
    if (n <= 2) return 0;
    
    // 只处理奇数：isPrime[i] 表示数字 (2*i+3) 是否为质数
    // 即 isPrime[0] -> 3, isPrime[1] -> 5, ...
    int size = (n - 1) / 2;
    boolean[] isPrime = new boolean[size];
    Arrays.fill(isPrime, true);
    
    int count = 1;  // 2 是质数
    
    for (int i = 0; i < size; i++) {
        if (isPrime[i]) {
            int p = 2 * i + 3;  // 对应的质数
            count++;
            
            // 从 p*p 开始标记奇数倍
            if ((long) p * p < n) {
                // p*p = (2i+3)^2 = 4i^2 + 12i + 9
                // 对应的索引为 (p*p - 3) / 2
                for (int j = (p * p - 3) / 2; j < size; j += p) {
                    isPrime[j] = false;
                }
            }
        }
    }
    
    return count;
}

线性筛（欧拉筛）​

埃氏筛会重复标记某些合数（例如 12 会被 2 和 3 分别标记一次）。线性筛保证每个合数只被标记一次，时间复杂度严格为 O(n)O(n)。

核心思想​

每个合数只被它的最小质因子标记。

/**
 * 线性筛（欧拉筛）
 * 时间复杂度: O(n)
 */
public int linearSieve(int n) {
    if (n <= 2) return 0;
    
    boolean[] isPrime = new boolean[n];
    Arrays.fill(isPrime, true);
    
    // 存储找到的质数
    List<Integer> primes = new ArrayList<>();
    
    for (int i = 2; i < n; i++) {
        if (isPrime[i]) {
            primes.add(i);
        }
        
        // 用已找到的质数标记合数
        for (int p : primes) {
            // 防止溢出
            if ((long) i * p >= n) break;
            
            isPrime[i * p] = false;
            
            // 关键：如果 p 整除 i，说明 i*p 已经被 p 标记过
            // 且 p 是 i*p 的最小质因子，停止标记
            if (i % p == 0) break;
        }
    }
    
    return primes.size();
}

为什么能保证 O(n)？​

关键在于 if (i % p == 0) break; 这行代码。当一个合数 x=i×px = i \times p 满足 p∣ip \mid i 时，pp 一定是 xx 的最小质因子。此时停止标记，后续更大的质数不会标记 xx。因此每个合数只被标记一次。

快速幂​

快速幂是一种高效计算 ana^n 的算法，将时间复杂度从 O(n)O(n) 降到 O(log⁡n)O(\log n)。

算法原理​

利用指数的二进制表示和平方的性质：

an=ab0⋅20⋅ab1⋅21⋅ab2⋅22⋯a^n = a^{b_0 \cdot 2^0} \cdot a^{b_1 \cdot 2^1} \cdot a^{b_2 \cdot 2^2} \cdots

其中 bib_i 是 nn 二进制表示的第 ii 位。

例如，a13=a11012=a8⋅a4⋅a1a^{13} = a^{1101_2} = a^8 \cdot a^4 \cdot a^1。

计算过程：

  1. 将指数 nn 写成二进制
  2. 从低位到高位遍历，如果当前位是 1，则乘上当前的 aa
  3. 每次迭代将 aa 平方，将 nn 右移一位

基本实现​

/**
 * 快速幂计算 a^n
 * 时间复杂度: O(log n)
 * 
 * @param a 底数
 * @param n 指数
 * @return a^n
 */
public long power(long a, int n) {
    long result = 1;
    
    while (n > 0) {
        // 如果当前位是 1，乘上当前的 a
        if ((n & 1) == 1) {
            result *= a;
        }
        // a 平方，n 右移
        a *= a;
        n >>= 1;
    }
    
    return result;
}

模幂运算​

在密码学和大数运算中，通常需要计算 an mod ma^n \bmod m。关键是在每次运算后取模，防止溢出。

/**
 * 快速幂取模：计算 (a^n) % mod
 * 时间复杂度: O(log n)
 * 
 * 应用场景：RSA 加密、费马小定理检验素数
 */
public long modPower(long a, long n, long mod) {
    long result = 1;
    a %= mod;  // 先取模，防止初始值过大
    
    while (n > 0) {
        if ((n & 1) == 1) {
            result = (result * a) % mod;
        }
        a = (a * a) % mod;
        n >>= 1;
    }
    
    return result;
}

矩阵快速幂​

矩阵快速幂可以高效计算斐波那契数列等递推关系：

/**
 * 矩阵快速幂计算斐波那契数列第 n 项
 * 时间复杂度: O(log n)
 */
public long matrixFibonacci(int n) {
    if (n <= 1) return n;
    
    // 斐波那契矩阵：[[1,1], [1,0]]
    long[][] matrix = {{1, 1}, {1, 0}};
    long[][] result = matrixPower(matrix, n - 1);
    
    return result[0][0];
}

/**
 * 矩阵快速幂
 */
private long[][] matrixPower(long[][] matrix, int n) {
    long[][] result = {{1, 0}, {0, 1}};  // 单位矩阵
    
    while (n > 0) {
        if ((n & 1) == 1) {
            result = multiply(result, matrix);
        }
        matrix = multiply(matrix, matrix);
        n >>= 1;
    }
    
    return result;
}

/**
 * 2x2 矩阵乘法
 */
private long[][] multiply(long[][] a, long[][] b) {
    long[][] c = new long[2][2];
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j];
        }
    }
    return c;
}

快速乘​

类似于快速幂的思想，可以优化乘法运算，常用于防止乘法溢出：

/**
 * 快速乘：计算 a * b % mod
 * 通过加法实现乘法，避免乘法溢出
 * 时间复杂度: O(log b)
 */
public long quickMultiply(long a, long b, long mod) {
    long result = 0;
    a %= mod;
    
    while (b > 0) {
        if ((b & 1) == 1) {
            result = (result + a) % mod;
        }
        a = (a + a) % mod;
        b >>= 1;
    }
    
    return result;
}

费马小定理与模逆元​

费马小定理​

如果 pp 是质数，且 gcd⁡(a,p)=1\gcd(a, p) = 1，则：

ap−1≡1(modp)a^{p-1} \equiv 1 \pmod{p}

由此可得：

ap−2≡a−1(modp)a^{p-2} \equiv a^{-1} \pmod{p}

这提供了一种使用快速幂求模逆元的方法。

模逆元计算​

/**
 * 费马小定理求模逆元
 * 前提：mod 是质数，且 gcd(a, mod) = 1
 */
public long fermatInverse(long a, long mod) {
    return modPower(a, mod - 2, mod);
}

/**
 * 扩展欧几里得求模逆元（通用方法）
 * 前提：gcd(a, mod) = 1
 */
public long extendedInverse(long a, long mod) {
    long[] result = extendedGcdLong(a, mod);
    return (result[1] % mod + mod) % mod;
}

private long[] extendedGcdLong(long a, long b) {
    if (b == 0) {
        return new long[]{a, 1, 0};
    }
    long[] temp = extendedGcdLong(b, a % b);
    long gcd = temp[0];
    long x = temp[2];
    long y = temp[1] - (a / b) * temp[2];
    return new long[]{gcd, x, y};
}

大数运算​

当数字超过 long 的范围时，需要使用字符串或数组来表示大数。

大数加法​

/**
 * 大数加法：字符串表示的大数相加
 */
public String addStrings(String num1, String num2) {
    StringBuilder result = new StringBuilder();
    
    int i = num1.length() - 1;
    int j = num2.length() - 1;
    int carry = 0;
    
    while (i >= 0 || j >= 0 || carry > 0) {
        int digit1 = i >= 0 ? num1.charAt(i--) - '0' : 0;
        int digit2 = j >= 0 ? num2.charAt(j--) - '0' : 0;
        
        int sum = digit1 + digit2 + carry;
        result.append(sum % 10);
        carry = sum / 10;
    }
    
    return result.reverse().toString();
}

大数乘法​

/**
 * 大数乘法：字符串表示的大数相乘
 * 时间复杂度: O(n * m)，其中 n 和 m 是两个数的长度
 */
public String multiplyStrings(String num1, String num2) {
    int n = num1.length();
    int m = num2.length();
    
    // 结果最多有 n + m 位
    int[] result = new int[n + m];
    
    // 从低位到高位计算
    for (int i = n - 1; i >= 0; i--) {
        for (int j = m - 1; j >= 0; j--) {
            int digit1 = num1.charAt(i) - '0';
            int digit2 = num2.charAt(j) - '0';
            
            int product = digit1 * digit2;
            int p1 = i + j, p2 = i + j + 1;  // 乘积的两个位置
            
            int sum = product + result[p2];
            result[p2] = sum % 10;
            result[p1] += sum / 10;
        }
    }
    
    // 构建结果字符串，跳过前导零
    StringBuilder sb = new StringBuilder();
    for (int digit : result) {
        if (!(sb.length() == 0 && digit == 0)) {
            sb.append(digit);
        }
    }
    
    return sb.length() == 0 ? "0" : sb.toString();
}

小结​

数学算法的核心要点：

算法时间复杂度核心思想
欧几里得算法O(log⁡min⁡(a,b))O(\log \min(a, b))gcd⁡(a,b)=gcd⁡(b,a mod b)\gcd(a, b) = \gcd(b, a \bmod b)
埃氏筛O(nlog⁡log⁡n)O(n \log \log n)标记质数的倍数
线性筛O(n)O(n)每个合数只被最小质因子标记一次
快速幂O(log⁡n)O(\log n)利用二进制分解指数

应用场景：

  • GCD/LCM：分数化简、周期计算
  • 质数筛：密码学、数论问题
  • 快速幂：RSA 加密、大数取模
  • 扩展欧几里得：解线性同余方程

练习​

基础：

  1. LeetCode 204：计数质数（埃氏筛）
  2. LeetCode 50：Pow(x, n)（快速幂）
  3. LeetCode 172：阶乘后的零（数学规律）
  4. LeetCode 191：位 1 的个数（位运算）

进阶： 5. LeetCode 372：超级次方（快速幂 + 递归） 6. LeetCode 365：水壶问题（GCD 应用） 7. LeetCode 149：直线上最多的点数（GCD 处理斜率） 8. LeetCode 326：3 的幂（数学方法）

挑战： 9. LeetCode 866：回文质数（筛法 + 构造回文） 10. LeetCode 908：最小差值 II（数学分析） 11. 实现 Miller-Rabin 素性测试 12. 实现大数阶乘计算