MillerRabin素数测试法

先给出wiki的数据

知道大家比较关心正确，网上教程好难找到数据，我就自己去wiki翻了
if n < 2,047, it is enough to test a = 2;
if n < 1,373,653, it is enough to test a = 2 and 3;
if n < 9,080,191, it is enough to test a = 31 and 73;
if n < 25,326,001, it is enough to test a = 2, 3, and 5
我自己验证了上面四组均通过

如果是uint范围：if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
如果是10^12范围：if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;
如果是10^14范围：if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
如果是10^18范围if n < 3,825,123,056,546,413,051, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, and 23.
如果是ull范围：if n < 18,446,744,073,709,551,616 = 2^64, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37.

给出一个MillerRabin的实现方式
我的思路大体是这样的，由费马小定理得 $a^{p-1}\equiv1 \pmod p$ap−1≡1(modp)，如果p是素数先要通过费马小定理的检测，然后再令$k=p-1$k=p−1，不断除2直到k为奇数，或者找到一个$a^k\equiv1 \pmod p$ak≡1(modp),但是$k/2$k/2不同余。
判断一下这个数是不是1或者p-1即可。

贴代码时间到

//请注意下这题的询问N是不正确的，而且还是一个ll范围的数。
//https://www.luogu.org/problemnew/show/U44392

#include<iostream>
#include<algorithm>
#include<ctime>
using namespace std;
#define ll long long

ll mymul(ll a, ll b, ll mod){
	ll ret = 0;
	while (b){
		if (b & 1) ret = (ret + a) % mod;
		a = (a + a) % mod;b >>= 1;
	}
	return ret;
}

ll mypow(ll a, ll b, ll mod) {
	ll ret = 1;
	while (b) {
		if (b & 1)ret = mymul(ret, a, mod);
		a = mymul(a, a, mod); b >>= 1;
	}
	return ret;
}
 
bool pd(ll a,ll p) {
	ll k = p - 1;
	if (mypow(a, k, p) != 1)return 0;
	while(k % 2 == 0 && mypow(a, k, p) == 1) k /= 2;
	ll tmp = mypow(a, k, p);
	return tmp == 1 || tmp == p - 1;
}

const int maxn = 100010;
bool vis[maxn];
int prime[maxn], cnt;
void init() {
	vis[1] = 1;
	for (int i = 2; i <maxn; i++) {
		if (!vis[i])prime[cnt++] = i;
		for (int j = 0; j < cnt&&prime[j] * i < maxn; j++) {
			vis[prime[j] * i] = 1;
			if (prime[j] % i == 0)break;
		}
	}
}
//一下依次是ull,10^18,10^14,10^12,uint范围
//ll pSet[12] = { 2,3,5,7,11,13,17,19,23,29,31,37 };
//ll pSet[9] = { 2,3,5,7,11,13,17,19,23 };
//ll pSet[7] = { 2,3,5,7,11,13,17 };
//ll pSet[4] = { 1,13,23,1662803 };
//ll pSet[3] = { 2,7,61 };

ll pSet[4] = { 2,13,23,1662803 };
const int num = 4;
bool rabin(ll p) {
	if (p < maxn) return !vis[p];
	for (int i = 0; i < num; i++) {
		if (!pd(pSet[i], p))return 0;
	}
	return 1;
}

int main() {
	init();
	//把最前面的垃圾数据读掉
	ll t; cin >> t;
	ll n;
	while (cin >> n) {
		if (rabin(n))cout << "Prime" << endl;
		else cout << "Not prime" << endl;
	}
}