How many prime numbers

题目链接

How many prime numbers

大素数的判定

解题思路

miller_rabin

miller_rabin 是一种概率性素数测试，原理基于费马素数测试，即如果 \(n\) 为素数，则在 \(1\sim n\) 中随机选择一个数 \(a\)，由费马小定理，一定有 \(a^{n-1}\equiv 1(\bmod n)\)，如果不满足则一定不是素数，但费马小定理的逆命题不一定成立，卡迈克尔数（费马伪素数）就是一种例外，但这样的数却并不多，miller_rabin 算法则是在费马素数检测的基础上尽可能能判断出卡迈克尔数，借助卡迈克尔数的一些性质，即卡迈克尔数一定由至少 \(3\) 个不同质数乘积组成，即卡迈克尔数一定不是 \(p^e\) 的形式，另外再由二次探测定理：\(\text { 如果 } p \text { 是奇素数，则 } x^2 \equiv 1(\bmod p) \text { 的解为 } x \equiv 1(\bmod p) \text { 或者 } x \equiv p-1(\bmod p)\)，其逆否命题为：如果 \(x^2 \equiv 1(\bmod p) \text { 的解都不为 } x \equiv 1(\bmod p) \text { 和 } x \equiv p-1(\bmod p)\)，则 \(p\) 不为奇素数，将 \(a^{n-1}\) 变换成 \(a^{u\times 2^t}=((a^{u})^{2^{t-1}})^2\)的形式，即计算 \(t-1\) 轮，判断 \((a^{u})^{2^{t-1}}\) 对 \(p\) 的模数是否为 \(1(\bmod p)\) 或 \(p-1(\bmod p-1)\)，如果都不是，说明其不是奇素数

• 时间复杂度：\(O(k\times logn)\)

代码

// Problem: hdu2138 How many prime numbers
// Contest: hdu
// http://acm.hdu.edu.cn/showproblem.php?pid=2138
// Memory Limit: 0 MB
// Time Limit: 1000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

// %%%Skyqwq
#include <bits/stdc++.h>

//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;

typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;

template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }

template <typename T> void inline read(T &x) {
    int f = 1; x = 0; char s = getchar();
    while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
    while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
    x *= f;
}

const int test_time=20;
int n;
LL x;
LL ksm(LL a,LL b,LL p)
{
	LL res=1%p;
	while(b)
	{
		if(b&1)res=(__int128)res*a%p;
		a=(__int128)a*a%p;
		b>>=1;
	}
	return res;
}
bool miller_rabin(LL n)
{
	if(n<3||n%2==0)return n==2;
	LL a=n-1,b=0;
	while(a%2==0)a>>=1,b++;
	for(int i=1,j;i<=test_time;i++)
	{
		LL x=rand()%(n-2)+2,v=ksm(x,a,n);
		if(v==1)continue;
		for(j=0;j<b;j++)
		{
			if(v==n-1)break;
			v=(__int128)v*v%n;
		}
		if(j>=b)return false;
	}
	return true;
}
int main()
{

    while(~scanf("%d",&n))
    {
        int res=0;
        while(n--)
        {
            scanf("%lld",&x);
            res+=miller_rabin(x);
        }
        printf("%d\n",res);
    }
    return 0;
}