每日算法 — 使用java实现Miller-Rabin:概率素性测试与密码学应用

引言

素性测试是数论与密码学的基础问题:给定一个大整数 n,如何高效判断它是素数还是合数?确定性的多项式时间算法(AKS)虽存在但效率较低,工程实践中广泛使用 Miller-Rabin 概率素性测试。该算法基于费马小定理与二次探测定理,几十轮测试即可将误判概率降至 2^(-80) 以下,被 OpenSSL、GnuPG 等主流密码库采用。

本文将用 Java 从零实现 Miller-Rabin 算法,从费马小定理、伪素数概念、二次探测定理到工程级实现细节,完整呈现这一经典算法的理论与工程实践。

数学基础

素数与合数

  • 素数:大于 1 且只有 1 和自身两个因子的正整数
  • 合数:大于 1 且不是素数的正整数

直接试除法时间复杂度为 O(√n),对于 1024 位大整数完全不可行。

费马小定理

若 p 是素数且 a 不是 p 的倍数,则:

a^(p-1) ≡ 1 (mod p)

逆否命题:若存在 a 使得 a^(n-1) ≢ 1 (mod n),则 n 是合数。

费马伪素数:满足费马小定理但实际是合数的 n 称为费马伪素数。

二次探测定理

若 p 是素数且 x² ≡ 1 (mod p),则 x ≡ ±1 (mod p)。

这是 Miller-Rabin 的核心:避免将卡迈克尔数(对所有 a 满足费马小定理的合数)误判为素数。

算法原理

Miller-Rabin 将 n-1 分解为 2^s × d(d 为奇数),然后对随机基数 a 进行测试:

  1. 计算 x = a^d mod n
  2. 若 x = 1 或 x = n-1,本轮通过
  3. 否则重复平方 x = x² mod n s-1 次
  4. 若任一次平方得到 n-1,本轮通过;否则 n 必为合数

正确性保证:合数 n 通过单轮测试的概率 ≤ 1/4;k 轮测试全部通过的概率 ≤ 4^(-k)。

基础实现

package com.algorithm.millerrabin;

import java.math.BigInteger;
import java.security.SecureRandom;

/**
 * Miller-Rabin 概率素性测试
 * 核心思路:将 n-1 分解为 2^s × d,对随机 a 测试
 *         a^d mod n 是否经过 n-1 或 1
 */
public class MillerRabin {
    private static final SecureRandom RANDOM = new SecureRandom();

    /**
     * 单轮 Miller-Rabin 测试
     * @return true 表示 n 通过本轮测试(可能为素数)
     */
    public static boolean millerRabinTest(BigInteger n, int s, BigInteger d, BigInteger a) {
        if (n.compareTo(BigInteger.ONE) <= 0) return false;
        if (n.equals(BigInteger.TWO) || n.equals(BigInteger.valueOf(3))) return true;
        if (n.mod(BigInteger.TWO).equals(BigInteger.ZERO)) return false;

        // 计算 a^d mod n
        BigInteger x = a.modPow(d, n);

        // 若 x = 1 或 x = n-1,通过
        if (x.equals(BigInteger.ONE) || x.equals(n.subtract(BigInteger.ONE))) {
            return true;
        }

        // 重复平方 s-1 次
        for (int r = 1; r < s; r++) {
            x = x.modPow(BigInteger.TWO, n);
            if (x.equals(n.subtract(BigInteger.ONE))) {
                return true;
            }
        }

        return false;  // n 是合数
    }

    /**
     * 多轮 Miller-Rabin 素性测试
     * @param n 待测试数
     * @param rounds 测试轮数(建议 ≥ 20)
     */
    public static boolean isPrime(BigInteger n, int rounds) {
        // 边界情况
        if (n.compareTo(BigInteger.TWO) < 0) return false;
        if (n.equals(BigInteger.TWO) || n.equals(BigInteger.valueOf(3))) return true;
        if (n.mod(BigInteger.TWO).equals(BigInteger.ZERO)) return false;

        // 分解 n - 1 = 2^s × d
        BigInteger nMinusOne = n.subtract(BigInteger.ONE);
        int s = nMinusOne.getLowestSetBit();  // 2^s 中最大的 s
        BigInteger d = nMinusOne.shiftRight(s);

        for (int i = 0; i < rounds; i++) {
            // 选随机基数 a ∈ [2, n-2]
            BigInteger a;
            do {
                a = new BigInteger(n.bitLength(), RANDOM);
            } while (a.compareTo(BigInteger.TWO) < 0 || 
                     a.compareTo(n.subtract(BigInteger.TWO)) > 0);

            if (!millerRabinTest(n, s, d, a)) {
                return false;  // 找到合数证据
            }
        }

        return true;  // 通过所有测试
    }

    /**
     * 默认轮数测试(40 轮,错误概率 < 2^-80)
     */
    public static boolean isPrime(BigInteger n) {
        return isPrime(n, 40);
    }
}

确定性 Miller-Rabin

对于 32 位以内的整数,可以使用一组固定的基数使 Miller-Rabin 成为确定性算法。这意味着测试这些范围内的任意整数,结果一定正确。

package com.algorithm.millerrabin;

import java.math.BigInteger;
import java.util.Arrays;
import java.util.HashSet;
import java.util.Set;

/**
 * 确定性 Miller-Rabin:对特定范围内的整数保证正确
 * 边界依据:Sorenson & Webster (2017) 证明 32 位以内可由前 13 个素数确定
 */
public class DeterministicMillerRabin {

    // 32 位以内的确定性基数集合
    private static final int[] DETERMINISTIC_WITNESSES_32 = {
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37
    };

    /**
     * 确定性素性测试(n < 3,317,044,064,679,887,385,961,981)
     * 即 2^82 以下的整数
     */
    public static boolean isPrime(BigInteger n) {
        if (n.compareTo(BigInteger.TWO) < 0) return false;
        if (n.equals(BigInteger.TWO) || n.equals(BigInteger.valueOf(3))) return true;
        if (n.mod(BigInteger.TWO).equals(BigInteger.ZERO)) return false;

        // 写 n - 1 = 2^s × d
        BigInteger nMinusOne = n.subtract(BigInteger.ONE);
        int s = nMinusOne.getLowestSetBit();
        BigInteger d = nMinusOne.shiftRight(s);

        for (int a : DETERMINISTIC_WITNESSES_32) {
            if (BigInteger.valueOf(a).compareTo(n) >= 0) break;

            BigInteger x = BigInteger.valueOf(a).modPow(d, n);
            if (x.equals(BigInteger.ONE) || x.equals(nMinusOne)) continue;

            boolean continueLoop = false;
            for (int r = 1; r < s; r++) {
                x = x.modPow(BigInteger.TWO, n);
                if (x.equals(nMinusOne)) {
                    continueLoop = true;
                    break;
                }
            }
            if (continueLoop) continue;

            return false;
        }

        return true;
    }
}

工程级优化

拒绝采样

先生成随机数,再判断其是否通过 Miller-Rabin。这样约一半的奇数会在 1 轮测试中被拒绝,极大减少计算量。

package com.algorithm.millerrabin;

import java.math.BigInteger;
import java.security.SecureRandom;

/**
 * 高效大素数生成器
 * 使用拒绝采样与轮数自适应
 */
public class PrimeGenerator {
    private static final SecureRandom RANDOM = new SecureRandom();

    /**
     * 生成指定位数的大素数
     */
    public static BigInteger generatePrime(int bits) {
        if (bits < 2) throw new IllegalArgumentException("位数必须 ≥ 2");

        // 自适应轮数:位数越多,所需轮数越多
        int rounds = computeRounds(bits);

        while (true) {
            // 生成随机候选数(最高位置 1 保证位数)
            BigInteger candidate = new BigInteger(bits, RANDOM);
            candidate = candidate.setBit(bits - 1);

            // 快速过滤:必须是奇数
            if (candidate.mod(BigInteger.TWO).equals(BigInteger.ZERO)) {
                continue;
            }

            // 试除小素数
            if (isDivisibleBySmallPrime(candidate)) {
                continue;
            }

            // Miller-Rabin 测试
            if (MillerRabin.isPrime(candidate, rounds)) {
                return candidate;
            }
        }
    }

    /**
     * 试除 2, 3, 5, 7, 11, ..., 997 等小素数
     * 90% 以上的合数会被快速过滤
     */
    private static final int[] SMALL_PRIMES = generateSmallPrimes(1000);

    private static boolean isDivisibleBySmallPrime(BigInteger n) {
        for (int p : SMALL_PRIMES) {
            if (n.mod(BigInteger.valueOf(p)).equals(BigInteger.ZERO)) {
                return n.equals(BigInteger.valueOf(p));
            }
        }
        return false;
    }

    private static int[] generateSmallPrimes(int limit) {
        boolean[] isPrime = new boolean[limit + 1];
        Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int i = 2; i * i <= limit; i++) {
            if (isPrime[i]) {
                for (int j = i * i; j <= limit; j += i) {
                    isPrime[j] = false;
                }
            }
        }
        java.util.List<Integer> primes = new java.util.ArrayList<>();
        for (int i = 2; i <= limit; i++) {
            if (isPrime[i]) primes.add(i);
        }
        return primes.stream().mapToInt(i -> i).toArray();
    }

    /**
     * 依据 NIST SP 800-22 推荐的轮数
     * 错误概率 < 2^-100
     */
    private static int computeRounds(int bits) {
        if (bits <= 512) return 7;
        if (bits <= 1024) return 15;
        if (bits <= 2048) return 20;
        if (bits <= 4096) return 25;
        return 40;
    }
}

强素数生成

RSA 等密码学应用需要”强素数”:p-1 与 p+1 都有大的素因子,避免 Pollard’s p-1 与 Williams’ p+1 攻击。

package com.algorithm.millerrabin;

import java.math.BigInteger;
import java.security.SecureRandom;

/**
 * 强素数生成:满足密码学额外安全要求
 */
public class StrongPrimeGenerator {
    private static final SecureRandom RANDOM = new SecureRandom();

    /**
     * 生成强素数 p,满足:
     * 1. p 是素数
     * 2. (p-1)/2 是素数(安全素数 / Sophie Germain 素数衍生)
     * 3. p-1 含有大的素因子
     * 4. p+1 含有大的素因子
     */
    public static BigInteger generateStrongPrime(int bits) {
        while (true) {
            BigInteger p = PrimeGenerator.generatePrime(bits);

            // 检查 (p-1)/2 也是素数
            BigInteger halfP = p.subtract(BigInteger.ONE).divide(BigInteger.TWO);
            if (!MillerRabin.isPrime(halfP, 20)) {
                continue;
            }

            // 检查 p-1 含有大素因子
            if (hasLargeFactor(p.subtract(BigInteger.ONE), bits / 2 - 10)) {
                return p;
            }
        }
    }

    /**
     * 检查 n 是否含有 ≥ minBits 位的素因子
     */
    private static boolean hasLargeFactor(BigInteger n, int minBits) {
        // 简化:检查 n 是否被某个大素数整除
        // 实际可使用 Pollard Rho 进行因数分解
        BigInteger minFactor = BigInteger.ONE.shiftLeft(minBits);
        for (int i = 0; i < 10; i++) {
            BigInteger a = new BigInteger(n.bitLength(), RANDOM).mod(n);
            BigInteger d = a.gcd(n);
            if (d.compareTo(BigInteger.ONE) > 0 && 
                d.compareTo(minFactor) >= 0) {
                return true;
            }
        }
        return false;
    }
}

误判概率分析

轮数 k 错误概率上界 备注
1 1/4 基础保证
5 1/1024 弱安全场景
20 < 10^(-12) 商业应用
40 < 10^(-24) 密码学推荐
80 < 2^(-160) NIST 最高级

理论上 Miller-Rabin 存在 Carmichael 数的反例(如 561 = 3×11×17 对所有互素基数都是费马伪素数),但通过二次探测可以正确识别它们。

性能基准

package com.algorithm.millerrabin;

import java.math.BigInteger;

/**
 * 性能基准测试
 */
public class PerformanceBenchmark {
    public static void main(String[] args) {
        int[] sizes = {256, 512, 1024, 2048};

        for (int bits : sizes) {
            long start = System.nanoTime();
            BigInteger prime = PrimeGenerator.generatePrime(bits);
            long cost = (System.nanoTime() - start) / 1_000_000;

            System.out.printf("生成 %d 位素数耗时: %d ms%n", bits, cost);
            System.out.printf("素数 (前 32 hex): %s...%n", 
                              prime.toString(16).substring(0, 32));
        }

        // 验证正确性
        BigInteger[] testCases = {
            BigInteger.valueOf(561),    // Carmichael 数(合数)
            BigInteger.valueOf(1729),   // Carmichael 数(合数)
            BigInteger.valueOf(1009),   // 素数
            BigInteger.valueOf(2).pow(31).subtract(BigInteger.ONE)  // Mersenne 素数
        };

        System.out.println("\n=== 正确性验证 ===");
        for (BigInteger n : testCases) {
            boolean result = MillerRabin.isPrime(n, 20);
            System.out.printf("%d: %s (判断: %s)%n", n, 
                              result ? "素数" : "合数", result);
        }
    }
}

复杂度分析

操作 时间复杂度 空间复杂度
单轮测试 O(k^3) O(k^2)
k 轮测试 O(k^4) O(k^2)
模幂运算 O(k^3) O(k^2)
大素数生成 O(k^4) 平均 O(k^2)

其中 k 为位数。

应用场景

  • RSA 密钥生成:需要大素数生成
  • Diffie-Hellman 密钥交换:素数模下的离散对数
  • ElGamal 签名:基于素数域
  • 椭圆曲线密码学:曲线参数选择
  • 伪随机数生成:素数序列
  • 哈希表大小选择:使用素数以减少冲突

与其他素性测试的对比

算法 类型 复杂度 实用性
试除法 确定性 O(√n) 仅小整数
Fermat 测试 概率 O(k^3) 易被卡迈克尔数欺骗
Solovay-Strassen 概率 O(k^3) 已较少使用
Miller-Rabin 概率 O(k^3) 工程首选
AKS 确定性 O(k^6) 理论重要
ECPP 概率 O(k^5) 速度极快

总结

Miller-Rabin 素性测试是数论与密码学结合的典范。短短几十行核心代码,即可在毫秒级判断数百位大整数是否为素数。掌握其原理与实现细节,不仅能帮助我们理解密码学底层机制,更能深入体会概率算法在工程中的精妙应用。

参考资料

  1. Miller, G. L. (1976). Riemann’s hypothesis and tests for primality. STOC ’75.
  2. Rabin, M. O. (1980). Probabilistic algorithm for testing primality. Journal of Number Theory.
  3. NIST SP 800-22: A Statistical Test Suite for Random and Pseudorandom Number Generators.
  4. Sorenson, J., & Webster, J. (2017). Strong pseudoprimes to twelve bases. Mathematics of Computation.