Miller-Rabin素性测试&Pollard-Rho

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/82314237

Code:

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
", max_fac);
	}
	return 0;
}

原文地址:https://www.cnblogs.com/JiaZP/p/14007712.html