【原创】质数专题筛质数 Miller_rabin Pollard_rho BZOJ3667 Rabin-Miller算法

本文链接：https://blog.csdn.net/c20182030/article/details/102653515

本文深入探讨了质数的定义、性质及筛选方法，详细介绍了埃拉托斯特尼筛法，并讨论了费马小定理与二次探测定理在质数判定中的应用。此外，文章还讲解了Pollard Rho算法在寻找合数因子中的高效运用。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

质数

说在前面

数论复习 Part 3。
我以前也写过一篇*儒雅随和*的质数专题，还故掉下篇，导致我始终不会米勒罗宾更不会Pollard Rho。
而这两年来我居然基本没见过哪里用过这两个算法。
这是因为你太菜了。

质数的定义和基本性质

不知道。

筛质数

太简单了，于是上代码。

const int MAXN=10050,MANA=10000;
bool vis[MAXN];
int cnt,pri[MAXN];

void shai()
{
	for(int i=2;i<=MANA;i++)
	{
		if(!vis[i]) pri[++cnt]=i;
		for(int j=1;j<=cnt && i*pri[j]<=MANA;j++)
		{
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}

（从我筛欧拉函数的代码里扣下来的）
我太简单了，所以我还要解释一波。

我们的代码流程是这样子的，对于每个数i，如果它未被筛走，则其为质数，存入质数表中。
然后把i与质数表里的每个数相乘，乘积必然为合数，就筛走它。
在枚举质数表里的质数的时候，从小到大枚举，如果i整除当前的质数，则不再往后枚举了。

这样就可以保证每个合数只会被它最小的质因子筛走了。
为什么？因为“如果i整除当前的质数（假设是j）”的话，那么i必然是j的若干次方，如果让i乘上接下来的质数（假设为k），得到的数的最小的质因子不是k而是i，因此就不要再往后了。
这样就可以保证每个合数只会被它最小的质因子筛走了。

质数判定

众所周知费马小定理， $p是质数\rArr \forall gcd(a,p)=1，a^{p-1}\%p=1 ~;~ \forall a\in N， a^p\%p=a$ ，如果我们找出很多a，来看 $a_i^{p}\%p$ 等不等于a，是不是就可以判断p是不是质数呢？

很遗憾， $p 是质数$ 只是 $\forall a\in N， a^{p}\%p=a$ 的充分条件，存在卡迈克尔数把它卡掉。
最小的卡迈克尔数是561=3*11*17，满足所有数的561次方模561等于其本身，与561互质的数的560次方模上561等于1。

所以我们再找一个定理，二次探测定理， $若p为质数,x^2\equiv 1\mod p,则x\equiv 1 \mod p或x\equiv p−1 \mod p$ 。

证明:
$x^2\equiv 1\mod p \iff x^2-1\equiv 0 \mod p \iff (x+1)(x-1) \equiv 0 \mod p \iff x\equiv 1 \mod p 或 x\equiv -1 \mod p$ 。

最后一步可能有些问题，万一 $p$ 不整除 $(x + 1)$ 也不整除 $(x - 1)$ ，但就是整除它们的积呢？
我们证明一下，设 $p\nmid(x+1)$ ，则有 $g c d (p, x + 1) = 1$ ， $\therefore 存在a,b使得ap+b(x+1)=1$ ，同时乘上 $(x - 1)$ 得 $a p (x - 1) + b (x + 1) (x - 1) = (x - 1)$ ，方程左边的两项都被 $p$ 整除，则 $p\mid右边$ 即 $p\mid (x-1)$ 。当 $p\nmid (x-1)$ 时同理。故 $p$ 总是整除 $(x + 1), (x - 1)$ 中的一个且仅有一个。
以上是我抄的欧几里得引理的证明流程，一个数整除两个互质的数的乘积则这个数整除这两个数中的一个。
（话说欧几里得在几何原本里面提出了多少引理定理怎么就这一个是欧几里得引理呢？）

那么二次探测原理对于我们的质数判定有什么用呢？

我们随机的找一个数a，如果 $a^{p-1}\mod p \neq 1$ ，那么 $p$ 肯定是合数，反之则有可能是质数。
然后要用二次探测原理是吧？ $x^2$ 哪儿来？我们知道，除了2以外的质数 $p$ 一定是奇数，这 $p - 1$ 不就是一个偶数吗？也就是说，我们已知 $({a^{\frac{p-1}{2}}})^{2}\equiv 1\mod p$ ，那么如果 $p$ 是质数的话，一定会有 $a^{\frac{p-1}{2}}\%p=\pm 1$ 。
如果等于 $1$ 的话，我们还可以试着再深入一些，如果 $\frac{p-1}{2}$ 还是偶数，我们就再用二次探测， $(a^{\frac{p-1}{4}})^2 \equiv 1\mod p\rArr a^{\frac{p-1}{4}}\% p=\pm 1$ 。以此类推。
如果 $\frac{p-1}{2}$ 是奇数，或者 $a^{\frac{p-1}{2}}\%p=-1$ ，则只能暂罢了，换下一个a。

另一种判定方法。
我们的上一种方法会跑几次？
我们把p-1分解为 $2^k\times d$ 的形式（当然 $d$ 是奇数），我们上一种方法跑 $k$ 次以后就会得到 $a^{2d}\% p=1\rArr a^d\%p=\pm1$ 这一步。
于是我们可以从 $2^d$ 开始，看它 $\%p$ 是不是等于 $\pm 1$ ，如果是则说明 $p$ 可能是质数，否则把 $2^d$ 平方，然后继续。
据说方法二会快一些。

至于选哪些a？
据说选 ${2,3,7,61,24251\}$ 的话 $10^{16}$ 以内只有 $46856248255981$ 会被误判为质数。
于是你可以多选几个， ${2,3,7,11,19,31,37,61\}$ ，或者我身边的dalao们都在用的东方的神秘数字 $19260817$ 。
（只有我一个人傻乎乎的写n以内的随机数？）

代码

#define ll long long
int base[10]={2,3,7,11,37,19260817};
bool Miller_Rabin(ll p)
{
    if(p<2) return 0;
	for(int i=0;i<=5;i++)
	{
		if(p==base[i]) return 1;
		if(p%base[i]==0) return 0;
	}
	ll d=p-1; int k=0;
	while((d&1)==0) d>>=1,k++;
	for(int z=0;z<=5;z++)
	{
		ll x=ksm(base[z],d,p),y;
		for(int i=1;i<=k && x!=1;i++,x=y)
			if((y=x*x%p)==1 && x!=p-1) return 0;
		if(x!=1) return 0;
	}
	return 1;
}

学习自firecrazy大佬，ksm是快速幂，其中x*x可能乘爆。
firecrazy说他之所以这么快，踩爆我们剩下的人，是神秘数字的功劳。

Pollard_Pho

这个算法能很快速的求出一个数的一个因数。
（1和这个数本身除外，平方根除外）
假设我们求的这个数是n。

怎么找？

方法一：for循环
方法二：随机生成数p，看n%p是否等于0。
这样的成功率是 $\frac{因数个数}{n}$ 的。

我们假设p是n的一个因数。
我们现在想要提高我们随机出p的概率。

生日悖论。

我们随机选出两个数a和b，它们的差值等于p的概率是多少？

举个具体的例子。

我们在1000个数里面随机选出k个数，其中某两个数之差（的绝对值）等于27的概率是多少？（27是我的数字desu）
k=2时，我们先选一个数i，那么下一个数就要选i+27或者i-27，这样的概率是 $\frac{2}{1000}$ 的。
k=3时，我们算不能找到的概率。我们先选了一个数i，就意味着有两个位置不能选了，下一步不选这两个位置的概率是 $\frac{998}{1000}$ ，此时又会有两个新的位置不能选，所以k=3的不能找到的概率是 $\frac{996*998}{1000^2}$ 。
综上，概率应为： $1-\frac{998*996*994*\cdots*(998-2*k+4)}{1000^{k-1}}$ 。
for循环一下发现，k=27的时候，成功率就大于一半了。k=67的时候，成功率就大于99%了。

其实我求错了，因为你可以选 $\bold{i+27,i-27}$ 而不选 $\bold i$ ，管他的，反正不是重点。

常见的生日悖论的描述方式是一个班里有生日相同的人的概率，这相当于把差值换成0，概率也更好算了。（就可以我的sb办法来做了）

总之结论是 $k=\sqrt{n}$ 的时候，成功率达到50%。

总之，我们就猜想，我们随机选k个数，然后其中某两个数的差的绝对值是n的因子。
这样当然可行，只是我们选 $k=\sqrt{n}$ 的时候，要 $k^2$ 的算差，时间复杂度还是 $\Theta(n)$ 的，白做了。

于是我们又猜想，我们随机选k个数，然后其中某个数的差值(设为 $p = ∣ a - b ∣$ )不与n互质。此时， $g c d (p, n)$ 就是n的一个因子。（废话，gcd能不是n的因子吗）
这样当然可行，只是我们要大概选 $\sqrt[4]{n}$ 个数，（我不知道为什么），这不优秀。

Pollard先生站出来说，我不要你随机，我要我随机，看我的随机数生成器：
$f(i)=(f(i-1)^2+p)\%n$
然后你选的k个数里面，第一个数任选，然后第二个数取 $f (1)$ ，……，第k个数取 $f (k - 1)$ 。

P先生，你的随机数生成器出锅了，你看，p取2，n=55的时候，你会出现循环。
啊？这是这个p的问题，你让他自己解决，就是说你换一个p试试看。
P先生，可是我怎么判断它是否出现循环啊。
啊？你连Floyd判环都不会吗？
P先生，我会Floyd求最短路……
啊？这是你的问题，你自己解决。

有环就会有什么？
就会有 $\colorbox{black}{\color{black}{情侣狗在夜深人静的时候走圈圈每天跑圈跑得气喘吁吁的时候你们这些闲庭信步的情侣狗就用看异类的眼神看着我虽然老子真的是异类虽然老子不稀罕但是啊啊啊}}$ 两个人竞争赛跑！（？）
竞争赛跑就会有什么？
就会有速度差！
有速度差就会有什么？
就会有超圈！
有速度差，有超圈就会有什么？
有环！

鼓掌！

所以我们每次取 $f (i)$ 和 $f (2 i)$ ，如果二者相等，就说明出现了超圈现象。

P先生啊，根据咕咕咕之家原理，你的随机数生成器是迟早要出锅的吧，特别是你的p是常数时候。
哦？感谢Pollard_rho算法的所有数组和函数们，我会反思自己生成随机数时出现的问题，希望所有人能越来越好。

还有大佬说什么路径倍增什么你跑到f(127)就可以停了。
我就完全不懂了，P先生也不懂。

以下的代码是BZOJ3667 Rabin-Miller算法的。

不打优化890ms。
打了floyd判环560ms。

代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;

#define ll long long
ll T,n,mx,base[10]={2,3,7,11,37,19260817};
ll abx(ll a){return a>=0?a:-a;}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll mul(ll a,ll b,ll c){return ((a*b-(ll)((long double)a/c*b+1e-8)*c)%c+c)%c;}
ll ksm(ll a,ll b,ll c)
{
	ll ans=1; a%=c;
	while(b)
	{
		if(b&1) ans=mul(a,ans,c);
		b>>=1,a=mul(a,a,c);
	}
	return ans;
}

bool Miller_Rabin(ll p)
{
    if(p<2) return 0;
	for(int i=0;i<=5;i++)
	{
		if(p==base[i]) return 1;
		if(p%base[i]==0) return 0;
	}
	ll d=p-1; int k=0;
	while((d&1)==0) d>>=1,k++;
	for(int z=0;z<=5;z++)
	{
		ll x=ksm(base[z],d,p),y;
		for(int i=1;i<=k && x!=1;i++,x=y)
			if((y=x*x%p)==1 && x!=p-1) return 0;
		if(x!=1) return 0;
	}
	return 1;
}

ll get_f(ll i,ll c,ll n){return (mul(i,i,n)+c)%n;}
ll Pollard_rho(ll n,ll c)
{
	ll one=get_f(0,c,n),two=get_f(one,c,n);
	while(one!=two)
	{
		ll d=gcd(n,abx(one-two));
		if(d!=1) return d;
		one=get_f(one,c,n),two=get_f(get_f(two,c,n),c,n);
	}
	return n;
}

void solve(ll n)
{
	if(n%2==0) {mx=max(mx,2ll); while(n%2==0) n/=2;}
	if(n==1 || Miller_Rabin(n)) {mx=max(mx,n); return;}
	ll tmp=n;
	while(tmp==n) tmp=Pollard_rho(n,rand()%(n-1)+1);
	while(n%tmp==0) n/=tmp;
	solve(tmp),solve(n);
}

int main()
{
	scanf("%lld",&T);
	while(T--)
	{
		scanf("%lld",&n);
		if(Miller_Rabin(n)) puts("Prime");
		else mx=0,solve(n),printf("%lld\n",mx);
	}
}