Miller-Rabin素性测试&Pollard-Rho
jiazhaopeng · · 个人记录
Miller-Rabin素性测试
分别用到费马小定理和二次探测定理
【算法流程】
（1）对于偶数和 0，1，2 可以直接判断。
（2）设要测试的数为 ，我们取一个较小的质数 ，设 ，满足 （其中 是奇数）。
（3）我们先算出 ，然后不断地平方并且进行二次探测（进行 次）。
（4）最后我们根据费马小定律，如果最后 ，则说明 为合数。
（5）多次取不同的 进行 素数测试，这样可以使正确性更高
————————————————
版权声明：本文为CSDN博主「forever_dreams」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/forever_dreams/article/details/82314237Code:
int pri[15] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
inline ll quickmul(ll x, ll y, ll M);
inline ll quickmi(ll x, ll k, ll M);
inline bool MR(ll x) {
if (x == 2) return true;
if (x % 2 == 0) return false;
if (x <= 1) return false;
//将x-1分解成(2^s)的形式
//随便找个素数a(<x)
//a^(x-1) -> a^((2^s)*t) -> (((a^t)^2)^2)... (mod x)
ll t = x - 1, s = 0;
while (!(t&1)) {
t >>= 1; s++;
}
for (register int i = 1; i <= 10 && pri[i] < x; ++i) {
int a = pri[i];
ll b = quickmi(a, t, x);//a^t
ll k;//tmp
for (register int j = 1; j <= s; ++j) {//s次平方
k = quickmul(b, b, x);
if (k == 1 && b != 1 && b != x - 1)//二次探测
return false;
b = k;
}
if (b != 1) return false;//费马小定理
}
return true;
}
int main() {
read(n);
if (MR(n)) puts("It is prime.");
else puts("It isn't prime.");
return 0;
}
Pollard-Rho
快速判断素数要用Miller-Rabin，而Pollard-Rho是用来分解质因数的。准确地说，一次PR是用来找随便一个因数的。复杂度玄学，大概是O(n^(1/4))。
这个博客看起来不错：Pollard-rho算法[因子分解算法]
对于模板题：P4718 【模板】Pollard-Rho算法，我们面对的其实不是纯粹的模板题，而是一个毒瘤卡常题，甚至我都不敢称之为卡常，它甚至要求我们去掉龟速乘的log，改为精度不高的long double快速乘。
面对这样的卡常题，我们需要具备必要的卡常技巧，对精度与时间复杂度的精确平衡，以及各种面向数据编程的能力。只可惜我在卡常方面还是太菜了，只会inline 和 register，因此以下代码我部分借鉴 lhm_ 大佬的代码。
由于为了卡常，该程序对Miller-Rabin做了修改，牺牲了很大的正确率，因此不卡常是还是建议写正常的Miller-Rabin写法。
inline ll quickmul(ll x, ll y, ll M) {
Re ll c = (long double)(x) * y / M + 0.5;
c = x * y - c * M;
return c < 0 ? c + M : c;
}
inline ll quickmi(ll x, ll k, ll M) {
Re ll res = 1;
while (k) {
if (k & 1) res = quickmul(res, x, M);
x = quickmul(x, x, M);
k >>= 1;
}
return res % M;
}
ll gcd(ll a, ll b) {
return b ? gcd(b, a % b) : a;
}
inline bool MR(ll x) {
if (x == 46856248255981ll || x == 1) return false;
if (x == 2 || x == 3 || x == 7 || x == 61 || x == 24251) return true;
Re ll b = 2;
Re ll k = x - 1;
while (k) {
Re ll cur = quickmi(b, k, x);
if (cur != 1 && cur != x - 1) return false;
if (k & 1 || cur == x - 1) break;
k >>= 1;
}
b = 61, k = x - 1;
while (k) {
Re ll cur = quickmi(b, k, x);
if (cur != 1 && cur != x - 1) return false;
if (k & 1 || cur == x - 1) return true;
k >>= 1;
}
return true;
}
ll max_fac;
ll F(ll num) {
if (num < 0) return -num;
return num;
}
ll f(ll x, ll c, ll M) {
return (quickmul(x, x, M) + c) % M;
}
inline ll PR(ll x) {
Re ll s = 0, t = 0, c = 1ll * rand() % (x - 1) + 1, val = 1;
for (register ll goal = 1; ; goal <<= 1, s = t, val = 1) {
for (register ll step = 1; step <= goal; ++step) {
t = f(t, c, x);
val = quickmul(val, F(t - s), x);
if (step % 127 == 0) {
Re ll d = gcd(val, x);
if (d > 1) return d;
}
}
Re ll d = gcd(val, x);
if (d > 1) return d;
}
return 1000000000;
}
void fac(ll x) {
if (x <= max_fac || x < 2) return ;
if (MR(x)) {
max_fac = x;
return ;
}
Re ll p = x;
while (p >= x) p = PR(x);
while (x % p == 0) x /= p;
fac(x); fac(p);
}
int main() {
read(t);
while (t--) {
read(n);
max_fac = 1;
fac(n);
if (max_fac == n) puts("Prime");
else printf("%lld\n", max_fac);
}
return 0;
}
常用模板
比较好在考场上写出来的板子：
inline ll Rand(ll mod) {
return ((1ull * rand() << 48) ^ (1ull * rand() << 32) ^ (1ull * rand() << 16) ^ rand()) % mod;//Bug
}
inline bool MR(ll P) {
if (P == 2 || P == 3 || P == 5 || P == 7) return true;
if (P <= 10) return false;
ll yu = P - 1;
int k = 0;
while (!(yu & 1)) yu >>= 1, ++k;
for (int _ = 1; _ <= 5; ++_) {
ll a = Rand(P - 1) + 1; a = quickpow(a, yu, P);
for (int i = 1; i <= k; ++i) {
ll y = Mul(a, a, P);
if (y == 1 && a != 1 && a != P-1) return false;//bug
a = y;
}
if (a != 1) return false;
}
return true;
}
ll jzpc, jzpP;
inline ll jzpRand(ll x) {
return (Mul(x, x, jzpP) + jzpc) % jzpP;
}
inline ll F(ll x) { return x < 0 ? -x : x; }
ll get_gcd(ll a, ll b) {
return b ? get_gcd(b, a % b) : a;
}
inline ll PR(ll n) {
jzpP = n; jzpc = Rand(n - 1) + 1;
for (ll s = 0, t = 0, goal = 1, mul = 1; ; goal <<= 1, s = t, mul = 1) {
for (int i = 1; i <= goal; ++i) {
t = jzpRand(t);
mul = Mul(mul, F(s - t), n);
if (i == 127) {
ll d = get_gcd(n, mul);
if (d > 1) return d;
}
}
ll d = get_gcd(n, mul);
if (d > 1) return d;
}
}
void Divi(ll n) {
if (MR(n)) return ++mp[n], void();
else {
ll d = PR(n);
while (d >= n) d = PR(n);
Divi(d), Divi(n / d);
}
}