【数学】Miller-Rabin算法素数测试

为了能够判断一个数是否是素数,我们很明显可以 O ( n ) 或 者 O ( n ) 打 表 O(\sqrt n)或者O(n)打表 O(n )O(n)但是在n巨大的时候这样太慢了,有没有更快的办法呢?

接下来,开始瞎搞。

如果要我在确定性算法和随机算法中做选择,当然要随机算法来瞎搞啦。
想一下,和质数有关的定理有什么?
首先我们知道,根据费马小定理,如果 p p p是一个质数那么对于任何a和p互质有 a p − 1 ≡ 1 ( m o d p ) a^{p-1}≡1(modp) ap11(modp)所以我们只要找到一个a,与那个p不满足上述条件,那么p就一定是合数,发过来呢?找不到p就是质数?不,反例很多,人们甚至给这样的反例取了个名字叫Carmichael数。
虽然Carmichael数并不多,但是仍然可以说明刚刚的方法好像并不靠谱,还有没有什么别的办法来升级呢?

我们很容易可以想到,对于方程 x 2 ≡ 1 ( m o d p ) x^2≡1(modp) x21(modp),首先这个方程应该有至少两个解 x ≡ ± 1 ( m o d p ) x≡±1(modp) x±1(modp),如果p真的是个质数,应该只有这两个解而已,否则可以知道
( x + 1 ) ( x − 1 ) ≡ 0 ( m o d p ) ⇒ p ∣ ( x + 1 ) 或 p ∣ ( x − 1 ) 或 ( p ∣ ( x + 1 ) 且 p ∣ ( x − 1 ) ) (x+1)(x-1)≡0(modp) \Rightarrow p|(x+1)或p|(x-1)或(p|(x+1)且p|(x-1)) (x+1)(x1)0(modp)p(x+1)p(x1)(p(x+1)p(x1))
如果 p ∣ ( x + 1 ) 且 p ∣ ( x − 1 ) p|(x+1)且p|(x-1) p(x+1)p(x1)
⇒ x + 1 = k p , x − 1 = j p ( k , j 均 为 常 数 ) \Rightarrow x+1=kp,x-1=jp(k,j均为常数 ) x+1=kp,x1=jp(k,j)
⇒ 2 = ( k − j ) p ( 直 接 由 两 式 相 减 得 ) \Rightarrow 2=(k-j)p(直接由两式相减得 ) 2=(kj)p()
如果 p p p是大于 2 2 2的质数,那么上面的式子显然不成立.
所以 p ∣ ( x + 1 ) 或 p ∣ ( x − 1 ) p|(x+1)或p|(x-1) p(x+1)p(x1)
最终得出 x ≡ ± 1 ( m o d p ) x≡±1(modp) x±1(modp)

因此如果方程 x 2 ≡ 1 ( m o d p ) x^2≡1(mod p) x21(modp)有一解 x 0 ∉ { 1 , − 1 } x_0\notin\{1,-1\} x0/{1,1}那么p不是素数。
因此,我们先看看之前的费马小定理 a p − 1 ≡ 1 ( m o d p ) a^{p-1}≡1(modp) ap11(modp)我们知道p-1是个偶数(2没必要这样判吧)所以设 p − 1 = 2 s ∗ d p-1=2^s*d p1=2sd那么
⇒ a 2 s ∗ d ≡ 1 ( m o d p ) \Rightarrow a^{2^s*d}≡1(modp) a2sd1(modp)
⇒ ( a 2 ( s − 1 ) ∗ d ) 2 ≡ 1 ( m o d p ) \Rightarrow (a^{2^{(s-1)}*d})^2≡1(modp) (a2(s1)d)21(modp)
诶!平方出来了,所以根据刚才的定理,我们知道 a 2 ( s − 1 ) ∗ d ≡ 1 ( m o d p ) 或 者 a 2 ( s − 1 ) ∗ d ≡ − 1 ( m o d p ) a^{2^{(s-1)}*d}≡1(modp)或者a^{2^{(s-1)}*d}≡-1(modp) a2(s1)d1(modp)a2(s1)d1(modp)如果 a 2 ( s − 1 ) ∗ d ≡ 1 ( m o d p ) a^{2^{(s-1)}*d}≡1(modp) a2(s1)d1(modp)那么我们将这个结果递归下去最终得到
⇒ a d ≡ 1 ( m o d p ) 或 a 2 r ∗ d ≡ − 1 ( m o d p ) \Rightarrow a^d≡1(modp)或a^{2^r*d}≡-1(modp) ad1(modp)a2rd1(modp)也就是说 a d a^d ad一开始就是1或者中间结果是-1
为了运算的方便,我们倒着来(毕竟求平方比开根方便)先把 a d a^d ad求出来再运算.
所以,如果n是质数那么 a d ≡ 1 ( m o d p ) 或 a 2 r ∗ d ≡ − 1 ( m o d p ) a^d≡1(modp)或a^{2^r*d}≡-1(modp) ad1(modp)a2rd1(modp)一定有一个成立.
由于可能是 a 2 r ∗ d ≡ − 1 ( m o d p ) a^{2^r*d}≡-1(modp) a2rd1(modp)的情况.
我们没有办法在一开始发现 a d ≡ 1 ( m o d p ) 或 a 2 r ∗ d ≡ − 1 ( m o d p ) a^d≡1(modp)或a^{2^r*d}≡-1(modp) ad1(modp)a2rd1(modp)不成立就判定这不是素数,需要继续运算。
下面是测一个a的情况的c++代码

bool R_M(long long d,long long a,long long M)//d什么的之前先处理了嘛
{
	long long x=ksm(a,d,M);
	if(x==1ll||x==M-1ll)//a^d≡1(modp)或a^{2^0*d}≡-1(modp)
		return 1;
	while(d<M-1)
	{
		x=ksc(x,x,M);//防止炸long long 的快速乘运算你可以理解为x=(x*x)%M;
		d<<=1ll;
		if(x==1ll)//这样就不可能之后有答案是-1了
			return 0;
		if(x==M-1ll)//找到了
			return 1;
	}
	return 0;//完了都没有-1一定不是质数
}

但是一个数还是有25%的概率骗过测试,所以我们只要多测几次就行了。
实际在竞赛中的应用,我们只需要检验一些固定的a即可。
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;
if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;
if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17;
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.
所以加上下面的代码,快速的质数判定就完成了

long long pr[20]={0,2,3,5,7,11,13,17,19,23,29,31};//随便整的,多几个没关系(如果卡常还是少几个吧)
bool Isprime(long long x)
{
	for(long long i=1;i<=MAXT;i++)
		if(pr[i]==x)//pr[i]是之前手打的几个质数,MAXT是测试组数
			return 1;
	if(x<=pr[MAXT])//我这里的质数是连续的才这样干,如果选的质数不连续不要这样干!
		return 0;
	long long d=x-1;
	while((d&1ll)==0)
		d>>=1ll;//求d
	for(long long i=1;i<=MAXT;i++)
		if(R_M(d,pr[i],x)==0)
			return 0;
	return 1;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值