Miller-Rabin素数测试

link

loj143

loj上板子题真难卡...

引入问题：给定一个数n，判断是不是质数。

这个问题很简单，可以在\(O(\sqrt n)\)内水过，不过如果毒瘤卡你时间，我们就需要更好的做法了。

Miller-Rabin就是很好的做法，可以在\(O(\log n)\)水过。

我们知道有个东西叫费马小定理，具体就是如果\(p\)是素数，那么有\(a^{p-1}\equiv 1\pmod p\)。它的逆命题，即如果\(a^{p-1}\equiv 1\pmod p\)，那么\(p\)是素数，显然我们随便举个例子就能说明这个命题fake。

我们还知道当\(p\)是素数且\(p>2\)的时候，对于方程\(x^2\equiv 1\pmod p\)，除了\(x\equiv1\pmod p\)或\(x\equiv-1\pmod p\)以外没有其它解。这个具体不太会证，好像跟二次剩余啥的有关。。

由于\(a^{p-1}\equiv 1\pmod p\)，我们把式子左右两边同时开平方，如果\(p\)是素数，肯定不会开出除了\(1\)和\(-1\)意外的其它数。

我们可以把\(p\)拆成\(2^kd+1\)的形式，我们先计算\(a^d\)，如果\(p\)是素数，那么在计算过程中，要么会出现\(-1\)，要么\(a^d\)一直是1。如果不满足这些条件，它肯定不是素数。

Miller-Rabin算法就是随机取一些\(a\)，都拿上面的方法判断即可。

loj数据太强悍，这里使用很多博客介绍的2,3,7,61,24251测试，还是WA了一半。所以我有加了10个[2,n-1]范围内的随机数才卡过去。

龟速乘太慢，这里偷懒写__int128

#include <cstdio>
#include <cstdlib>
#include <ctime>
using namespace std;

long long turtle(long long x, long long y, long long p)
{
	return ((__int128)x * (__int128)y % p);
	// long long ans = 0;
	// for (x %= p; y > 0; x = (x + x) % p, y >>= 1) if (y & 1) ans = (ans + x) % p;
	// return ans;
}

long long qpow(long long x, long long y, long long p)
{
	long long ans = 1;
	for (x %= p; y > 0; x = turtle(x, x, p), y >>= 1) if (y & 1) ans = turtle(ans, x, p);
	return ans;
}

bool Miller_Rabin(long long n, long long a)
{
	if (n == a) return true;
	if (n < 2 || n % a == 0) return false;
	long long d = n - 1;
	while ((d & 1) == 0) d >>= 1;
	long long t = qpow(a, d, n);
	while (d != n - 1 && t != 1 && t != n - 1) t = turtle(t, t, n), d <<= 1;
	return (t == n - 1) || ((d & 1) == 1); //对于计算到-1即可停下返回真，或一开始就是1
	//如果没有经过-1而直接计算到了1，或者全程没有经过-1，说明不满足性质，一定是合数
}

bool zz(long long x)
{
	if (x == 2) return true;
	for (int i = 1; i <= 10; i++)
		if (Miller_Rabin(x, rand() % (x - 2) + 2) == false) return false;
	return true;
}

bool prime(long long x)
{
	return Miller_Rabin(x, 2) && Miller_Rabin(x, 3) && Miller_Rabin(x, 7) && Miller_Rabin(x, 61) && Miller_Rabin(x, 24251) && zz(x);
}

int main()
{
	srand(time(0));
	long long x;
	while (scanf("%lld", &x) == 1)
	{
		if (prime(x)) printf("Y\n");
		else printf("N\n");
	}
	return 0;
}