Miller-Rabin素数测试算法

发表于: 2025-05-16 更新于: 2025-05-16 浏览: 1

title: Miller-Rabin素数测试算法
date: 2022-03-15 18:42:01
categories: 算法
tags: [素数,算法]

mathjax: true

学习RSA过程中用到了大素数，记录学习过程。

目前没有算法能够直接生成一个素数（据本人了解），只能通过不断生成随机数直到产生素数的方法来生成。如果$n$较小，可以直接对$2$到$\sqrt{n}$ 的所有整数取模验证，但当$n$特别大时，$O(\sqrt{n})$ 复杂度耗费的大量时间是不能被接受的。此时如果牺牲部分准确度来换取时间的优化，即Miller-Rabin算法的思想。

根据百度百科，使用快速幂运算，这个算法的时间复杂度是 $O(𝑘\log^{3𝑛)$ ，判错概率是$4}{-k}$，$k$为判断次数。下面为该算法的实现原理。

引理

1. 费马小定理

如果$p$是素数，且 $gcd(a, p) = 1$ ，则有 $a^{p-1}\equiv 1 \pmod {p}$ 。

它的逆命题并不能推出$p$为素数，但如果不满足 $a^{p-1} \equiv 1 \pmod {p}$ 的$p$一定是合数。

2. 有限域上的平方根定理

如果$p$是一个奇质数，且 $e \geq 1$ ，则方程

$$
x^{2 \equiv 1 \pmod {p}e}
$$

仅有两个根 $x=1$ 或者 $x=-1$ ，$\pm1$也称为1的平凡平方根。

在 $e=1$ 即模$p$的情况下，$x=-1$等价于$x=p-1$ 。

把 $x = 1$ 和 $x = p-1$ 称为“$x$对模$p$来说1的平凡平方根”不太通顺，Miller-Rabin索性测试用到这个方程：$x^{2 \equiv 1 \pmod {n}$ 。如果一个数$x$满足方程$x}2 \equiv 1 \pmod {𝑛}$，但$x$不等于平凡平方根1或$n-1$ ，那么称$x$是对模n来说1的“非平凡”平方根。例如，$x=6$ , $n=35$ ，6是对模35来说1的非平凡平方根。

推论：如果对模n存在1的非平凡平方根，则n是合数。

Miller-Rabin算法

判断一个大数n是否为素数时，可以利用上面的费马小定理 $a^{n-1} \equiv 1 \pmod {n}$ ，如果不满足该式，则可断定n为合数。但在该判断前可以用引理2的推理进一步增强准确性。

对于偶数$n-1$ ，一定可以拆分为 $n-1 = 2 ^{s * d$ ：

$$

a} {n-1} = a ^{2 {s} * d} = (((a^d)2)\cdots)^{2

$$

从$x=a}d$开始，依次平方s次，每次平方的时候模n。由有限域上的平方根定理可知，如果模n时结果为1，那么$x$一定是1或者$n-1$，如果不满足则n为合数。如果$x$为1的平凡平方根或模n的结果不为1，$x = x^2$继续下一次平方，然后模n判断。

当平方到了s次，此时$x = a ^ {n-1}$ ，判断$x$模n的结果是否为1。如果不为1判定n为合数，如果为1则可认为n大概率为素数。

算法的代码实现

在进行$a^{{d} \pmod n$ 和 $x}2 \pmod n$时如果直接运算，大概率会溢出，所以需要使用快速幂和快速乘。

#include <stdio.h>
#include <stdlib.h>

typedef long long ll;

ll qMul(ll a, ll b, ll mod);

ll qPow(ll base, ll power, ll mod);

int Miller_Rabin(ll n, int repeat);

int main() {
    ll n;
    while (printf("Input a n to judge:\n") && scanf("%lld", &n) != EOF) {
        printf(Miller_Rabin(n, 10) ? "Yes\n\n" : "No\n\n");
    }
    return 0;
}

//quick multiply
ll qMul(ll a, ll b, ll mod) {
    ll ret = 0;
    while (b) {
        if (b & 1) ret = (ret + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return ret;
}

//quick power
ll qPow(ll base, ll power, ll mod) {
    ll ret = 1;
    while (power) {
        if (power & 1) ret = qMul(ret, base, mod);
        base = qMul(base, base, mod);
        power >>= 1;
    }
    return ret % mod;
}

int Miller_Rabin(ll n, int repeat) {
    if (n == 2 || n == 3) return 1;
    if (n % 2 == 0 || n == 1) return 0;

    ll d = n - 1;
    int s = 0;
    while (!(d & 1)) ++s, d >>= 1;

    for (int i = 0; i < repeat; i++) {
        ll a = rand() % (n - 3) + 2;
        ll x = qPow(a, d, n);
        ll y = 0;
        for (int j = 0; j < s; j++) {
            y = qMul(x, x, n);
            if (y == 1 && x != 1 && x != (n - 1)) return 0;
            x = y;
        }
        if (y != 1) return 0;
    }

    return 1;
}