Fermat素性检验

费马小定理：$n$是质数，$gcd\left(a,n\right)=1$,则$a^{n-1}\equiv 1\left(\mathrm{mod}n\right)$
反过来不一定成立，即：如果$a^{n-1}\not\equiv 1\left(\mathrm{mod}n\right)$,则$n$也不一定是合数，例如$2^{340}\equiv 1\left(\mathrm{mod}341\right)$

所以这里的思想就是，在$\left[2,n-1\right]$中随机挑选$a$，检测是否每一次都有$a^{n-1}\equiv 1\left(\mathrm{mod}n\right)$
如果不是，则$n$是合数（因为此时$\forall a\in \left[2,n-1\right],s.t.gcd\left(a,n\right)=1$, 如果$n$是质数，则由费马小定理，矛盾）

#include <cstdio>
#include <stdlib.h>
typedef long long LL;
LL quick_pow(LL a, LL b, LL p) {
	LL res = 1;
	while (b) {
		if (b & 1)res = (res * a) % p;
		a = (a * a) % p;
		b >>= 1;
	}
	return res;
}
const int test_time = 500;
bool judge(int n) {
	if (n < 3)return n == 2;
	for (int i = 1; i <= test_time; ++i) {
		int a = rand() % (n - 2) + 2;
		if (quick_pow(a, n - 1, n) != 1)return false;
	}
	return true;
}

Carmichael number

$n$是一个正合数，且$\forall b,gcd\left(b,n\right)=1$,有$b^{n-1}\equiv 1\left(\mathrm{mod}n\right)$，则称$n$是一个卡迈克尔数/费马伪素数/Carmichael number

性质

性质1：Carmichael number必定是至少三个素数的乘积

Korselt’s criterion

$n$是一个正合数，当且仅当对于所有的素因数$p$，$p^2\nmid n$且$(p-1)\mid (n-1)$
其实这也说明了Carmichael number必然是一个奇数

例题

https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&category=12&page=show_problem&problem=947

#include<cstdio>
const int N = 65000;
bool visit[N + 5];
int primes[N / 3], cnt;
void init(int n) {
	primes[0] = 2;
	cnt = 1;
	for (int i = 4; i <= n; i += 2) {
		visit[i] = true;
	}
	for (int i = 3; 3 * i <= n; i += 2) {
		if (!visit[i]) {
			primes[cnt] = i;
			++cnt;
		}
		for (int j = i * i; j <= n; j += i) {
			visit[j] = true;
		}
	}
}
bool judge(int n) {
	//n是一个奇合数
	if ((n & 1) == 0 || !visit[n])return false;
	// 遍历n的素因子
	// n是一个奇数，所以最大的因子不超过n/3
	for (int i = 1; primes[i] * 3 <= n; ++i) {
		if (n % primes[i] != 0) {
			continue;
		}
		//p^2 \not\mid n and (p-1)|(n-1)
		if (n % (primes[i] * primes[i]) == 0 || (n - 1) % (primes[i] - 1) != 0) {
			return false;
		}
	}
	return true;
}
int main() {
	init(N);
	int n;
	while (scanf("%d", &n), n) {
		if (judge(n))printf("The number %d is a Carmichael number.\n", n);
		else printf("%d is normal.\n", n);
	}
	return 0;
}

Miller-Rabin素性测试

引理

如果$p$是质数，则$x^2\equiv 1\left(\mathrm{mod}p\right)$的解为$x\equiv 1\left(\mathrm{mod}p\right)$或$x\equiv -1\left(\mathrm{mod}p\right)$
证明：$p\mid \left(x^2-1\right)\Rightarrow p\mid \left(x-1\right)\left(x+1\right)$
又因为$p$是质数，$p\mid \left(x-1\right)$或$p\mid \left(x+1\right)$

原理

结合引理和Fermat素性检验

$n$是奇质数
考虑$a<n$,
$n-1$是一个偶数，设$n-1=2^{s}u$
如果$a^u\equiv 1\left(\mathrm{mod}n\right)$,或者$\exists t<s,$满足$a^{2^{t}u}\equiv -1\left(\mathrm{mod}n\right)$ ，则$n$可能是质数，否则就是合数

证明：
当$a^u\equiv \pm 1\left(\mathrm{mod}n\right)$,则$a^{n-1}={\left(a^u\right)}^{2^s}\equiv 1\left(\mathrm{mod}n\right),则$n$有可能是质数，返回

当$a^u\not\equiv \pm 1\left(\mathrm{mod}p\right)$，我们接着平方
如果$a^{2u}\equiv 1\left(\mathrm{mod}n\right)$，因为$a^u\not\equiv \pm 1\left(\mathrm{mod}p\right)$，根据引理，$n$是合数
如果$a^{2u}\equiv -1\left(\mathrm{mod}n\right)$,如果$2u=n-1$,则根据费马小定理$n$是合数，否则$n$有可能是质数，返回
如果$a^{2u}\not\equiv \pm 1\left(\mathrm{mod}n\right)$，我们就接着平方

同理对于$t<s$
如果$a^{2^{t}u}\equiv 1\left(\mathrm{mod}n\right)$，因为$a^{2^{t-1}u}\not\equiv \pm 1\left(\mathrm{mod}p\right)$，根据引理，$n$是合数
如果$a^{2^{t}u}\equiv -1\left(\mathrm{mod}n\right)$,$n$有可能是质数，返回
如果$a^{2^{t}u}\not\equiv \pm 1\left(\mathrm{mod}n\right)$，我们就接着平方

到最后$a^{n-1}\equiv 1\left(\mathrm{mod}n\right)$则$n$是合数

此外
针对$\left[1,2^{32}\right)$范围内的数，只需要选择 { 2 , 7 , 61 } \left\{2,7,61\right\} {2,7,61}作为基底进行Miller-Rabin素性测试就可以确定素性
针对$\left[1,2^{64}\right)$范围内的数，只需要选择 { 2 , 325 , 9375 , 28178 , 450775 , 9780504 , 1795265022 } \left\{2, 325, 9375, 28178, 450775, 9780504, 1795265022\right\} {2,325,9375,28178,450775,9780504,1795265022}作为基底进行Miller-Rabin素性测试就可以确定素性
针对$\left[1,2^{78}\right)$范围内的数，只需要选择 { 2 , 3 , 5 , 7 , 11 , 13 , 17 , 19 , 23 , 29 , 31 , 37 } \left\{2,3,5,7,11,13,17,19,23,29,31,37\right\} {2,3,5,7,11,13,17,19,23,29,31,37}(前12个质数)作为基底进行Miller-Rabin素性测试就可以确定素性
注意：
如果用这种方法，
1.所有数都要取一遍，不能只选择小于$n$的
2.$a$换成$a\mathrm{mod}n$
3.如果$a\equiv 0\left(\mathrm{mod}n\right)$则直接通过该轮测试

另外，假设 广义 Riemann 猜想(generalized Riemann hypothesis, GRH)成立，则对数$n$最多只需要测试$\left[2,\min \left\{n-2,\lfloor 2{\ln }^{2}n\rfloor \right\}\right]$中的全部整数即可确定$n$的素性

#include <cstdio>
#include <stdlib.h>
typedef long long LL;
LL quick_pow(LL a, LL b, LL p) {
	LL res = 1;
	while (b) {
		if (b & 1)res = (res * a) % p;
		a = (a * a) % p;
		b >>= 1;
	}
	return res;
}
const int test_time = 500;
bool Rabin_Miller(int n) {
	if (n < 3 || (n & 1) == 0)return n == 2;
	int u = n - 1, t = 0;
	// n - 1 = u * (2^t)
	while ((u & 1) == 0) {
		u >>= 1;
		++t;
	}
	for (int i = 1; i <= test_time; ++i) {
		int a = rand() % (n - 2) + 2, v = quick_pow(a, u, n);
		if (v == 1)continue;
		int s;
		for (s = 0; s < t; ++s) {
			//a^{u * 2^s}= -1 (mod n)
			if (v == n - 1)break;
			v = 1LL * v * v % n;
		}
		if (s == t)return false;
	}
	return true;
}

第二种：$\left[1,2^{32}\right)$

#include <cstdio>
#include <cstdlib>
#include <vector>
using namespace std;
typedef long long LL;
LL quick_pow(LL a, LL b, LL p) {
	LL res = 1;
	while (b) {
		if (b & 1)res = (res * a) % p;
		a = (a * a) % p;
		b >>= 1;
	}
	return res;
}
vector<unsigned> bases1 = { 2,7,61 };
bool Rabin_Miller2(unsigned n) {
	if (n < 3 || (n & 1) == 0)return n == 2;

	unsigned u = n - 1, t = 0;
	// n - 1 = u * (2^t)
	while ((u & 1) == 0) {
		u >>= 1;
		++t;
	}
	for (unsigned a : bases1) {
		a %= n;
		if (a == 0)continue;
		unsigned v = quick_pow(a, u, n);
		if (v == 1)continue;
		unsigned s;
		for (s = 0; s < t; ++s) {
			//a^{u * 2^s}= -1 (mod n)
			if (v == n - 1)break;
			v = 1ULL * v * v % n;
		}
		if (s == t)return false;
	}
	return true;
}

第三种：$\left[1,2^{64}\right)$
__int128要gcc或者g++才有，如果你能确保不溢出，换成long long也可以

#include <cstdio>
using namespace std;
typedef long long LL;
LL quick_pow(LL a, LL b, LL p) {
	LL res = 1;
	while (b) {
		if (b & 1)res = (res * a) % p;
		a = (a * a) % p;
		b >>= 1;
	}
	return res;
}
const LL bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
bool Rabin_Miller(LL n) {
    if (n < 3 || (n & 1) == 0)return n == 2;

    LL u = n - 1, t = 0;
    // n - 1 = u * (2^t)
    while ((u & 1) == 0) {
        u >>= 1;
        ++t;
    }
    for (int i = 0; i < 7; ++i) {
        LL a = bases[i] % n;
        if (a == 0)continue;
        LL v = quick_pow(a, u, n);
        if (v == 1)continue;
        LL s;
        for (s = 0; s < t; ++s) {
            //a^{u * 2^s}= -1 (mod n)
            if (v == n - 1)break;
            v = (__int128)v * v % n;
        }
        if (s == t)return false;
    }
    return true;
}

第四种：$\left[1,2^{78}\right)$
__int128要gcc或者g++才有，如果你能确保不溢出，换成long long也可以

#include <cstdio>
using namespace std;
typedef long long LL;
LL quick_pow(LL a, LL b, LL p) {
	LL res = 1;
	while (b) {
		if (b & 1)res = (res * a) % p;
		a = (a * a) % p;
		b >>= 1;
	}
	return res;
}
const LL bases[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
bool Rabin_Miller(LL n) {
    if (n < 3 || (n & 1) == 0)return n == 2;

    LL u = n - 1, t = 0;
    // n - 1 = u * (2^t)
    while ((u & 1) == 0) {
        u >>= 1;
        ++t;
    }
    for (int i = 0; i < 12; ++i) {
        if(bases[i] == n) return true;
        LL a = bases[i] % n;
        if (a == 0)return false;
        LL v = quick_pow(a, u, n);
        if (v == 1)continue;
        LL s;
        for (s = 0; s < t; ++s) {
            //a^{u * 2^s}= -1 (mod n)
            if (v == n - 1)break;
            v = (__int128)v * v % n;
        }
        if (s == t)return false;
    }
    return true;
}

参考：
https://www.ams.org/journals/mcom/1990-55-191/S0025-5718-1990-1023756-8/
https://miller-rabin.appspot.com/#