Miller-Rabin 素数判定 与 较大数取模

米勒拉宾判定

LOJ – 质数判定
最近学了这个随机化的素数判定方法,终于学懂了(我之前学的肯定是假的
给一个建议:最好不要去Wikipedia – 米勒-拉宾素性检验上看
我随便讲一下,因为最近有点事多
首先对于一个数 $p$ ,如果它是素数,那么由于费马小定理:
$$ a^{p-1}\equiv 1\pmod{p}$$
看上去似乎只要枚举一些 $a$ 判断是否满足就行了,但实际上很容易被卡掉
这时候考虑引入其他的限制来进行进一步的判断

若 $p$ 是奇素数,则 $a^2\equiv 1\pmod{p}$ 的解只可能为 $a\equiv 1$ 或 $a\equiv -1$
那么考虑将偶数 $n=p-1$ 拆分成 $n=2^r\times d$ ,那么若存在 $1\leq s\leq r$ ,有
$$
\begin{eqnarray}
a^{2^s d}&\equiv& 1\\
a^{2^{s-1} \ \ d}&\not \equiv& -1, 1
\end{eqnarray}
$$
则 $p$ 不是一个质数
当然,若存在 $a$ ,使
$$ a^{p-1}\not \equiv 1\pmod{p} $$
那么 $p$ 同样不是一个质数
正确率好像很高,但是我懒得写了
复杂度 $O(\log{p})$

较大数取模

发现 $p$ 可以到 $10^{18}$ ,那么可能要用慢速乘,但是这样一次乘法要 $\log{p}$ ,会跑得很慢,所以要寻找一个较大数取模的快速方法
由于
$$a\times b\mod p=a\times b-\lfloor\frac{a\times b}{p}\rfloor\times p$$
使用 long double 来计算 $\lfloor\frac{a\times b}{p}\rfloor$ ,误差较小, long double 存不下就会舍弃低位,再转成 long long ,小于 $2^{63}$ 的位不会挂,而模数小于 $2^{63}$ ,于是可以快速的求出取模后的数
大概实现方法

inline ll mul (ll x, ll y, ll MOD) {
    ll tmp = x * y - (ll)((long double)x / MOD * y + 0.5) * MOD;
    return tmp < 0 ? tmp + MOD : tmp;
}

程序实现

给个板子

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>

#define reg register
typedef long long ll;

ll N;
ll primes[9] = {3, 5, 7, 11, 13, 17, 19, 23};

inline ll mul (ll x, ll y, ll MOD) {
    ll tmp = x * y - (ll)((long double)x / MOD * y + 0.5) * MOD;
    return tmp < 0 ? tmp + MOD : tmp;
}

inline ll quick_pow (ll x, ll p, ll MOD) {
    ll ans = 1;
    while (p) {
        if (p & 1) ans = mul(ans, x, MOD);
        x = mul(x, x, MOD);
        p >>= 1;
    }
    return ans;
}

inline bool isPrime (ll P) {
    if (P < 2) return false;
    if (P == 2) return true;
    if (P & 1 == 0) return false;
    ll N = P - 1;
    int t = 0;
    while (!(N & 1LL)) N >>= 1LL, t++;
    for (int T = 0; T < 8; ++T) {
        if (P == primes[T]) return true;
        ll x = quick_pow(primes[T], N, P), y = x;
        for (int i = 0; i < t; ++i) {
            x = mul(x, x, P);
            if (x == 1 && !(y == 1 || y == P - 1)) return false;
            y = x;
        }
        if (x != 1) return false;
    }
    return true;
}

int main () {
    while (read(N)) {
        if (isPrime(N))
            puts("Y");
        else
            puts("N");
    }
    return 0;
}

引用资料

sys_con

OIer,常规 / 竞赛都渣得不行

发表评论

电子邮件地址不会被公开。 必填项已用*标注