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