【数论】Miller_Rabin算法

最新推荐文章于 2022-04-19 19:27:23 发布

lazy-sheep

最新推荐文章于 2022-04-19 19:27:23 发布

阅读量192

点赞数 2

分类专栏：算法数论文章标签：数论

本文链接：https://blog.csdn.net/zsjzliziyang/article/details/99763123

版权

算法同时被 2 个专栏收录

24 篇文章 2 订阅

订阅专栏

数论

12 篇文章 0 订阅

订阅专栏

由费马小定理：

若p是质数，p不整除于a，则

$a^{p-1}\equiv 1(mod~p)$

猜想：

若p不整除于a,且 $a^{p-1}\equiv 1(mod~p)$ ，则p是质数

显然不一定成立，但多试几个a，好像还挺准的？

于是我们可以得出一个伪Miller_Rabin算法

对于a={2、3、5、7、11、13、17、19、23、31、101}

若都满足 a=p 或 $a^{p-1}\equiv 1(mod~p)$ ，且 $p\neq 1$ 则p是质数。

实测证明，这种算法在10^18以内不会被卡。

但是严谨的我们并不能满足于此，如何优化？

二次探测定理：

若 $a^2\equiv 1(mod~p)&&(a\neq 1,a\neq p -1)$ 时，p一定是合数。

在p是质数的情况下，一定有 $a^2\equiv 1(mod~p)$ 的解为a=1或a=p-1

证明：

$a^2-1\equiv 0~(mod~p)$

$(a+1)(a-1) \equiv 0~(mod~p)$ $a^2\equiv 1~(mod~p)$

解得：

a=1~||~a=-1

即

a=1~||~a=p-1

则若 $a\neq 1\cup a\neq p -1$ ，p一定不是质数

我们可以在计算 $a^{p-1}$ 时验证这条式子

于是我们就得出了专业的Miller_Rabin算法

a取若干个素数，我们要判断 a=p 或 $a^{p-1}\equiv 1(mod~p)$ ，若所有的a都通过测试，那么则称p是质数。

设 p-1=2^kd ，我们先计算 a^d ，再讲得到的值平方k次。

若 a^d =1则直接通过测试（也可以不用判断）

接下来在平方的过程中，若如果平方出来的数=1且被平方的数不等于1和p-1那么直接退出，p不是质数。

最后再判断一次 $a^{p-1}\equiv 1(mod~p)$ ，若不满足也不是质数。

例题：

4278. 【NOIP2015模拟10.29B组】质数

Code

#include<cstdio> 
#include<algorithm>
#include<cstring>
#include<cstdlib>
#define I int
#define F(i,a,b) for(I i=a;i<=b;i++)
#define Fd(i,a,b) for(I i=a;i>=b;i--)
#define mem(a,b) memset(a,b,sizeof(a))
#define N 10000000
#define ll long long
#define P(x) freopen(x".in","r",stdin);freopen(x".out","w",stdout)
#define rt return
using namespace std;
void rd(ll &x){
	x=0;char ch=getchar();
	while(ch<'0'||ch>'9') ch=getchar();
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
}
ll x,n;
I a[13]={2,3,5,7,11,13,17,19,23,29,31,37,101};
ll ksc(ll x,ll y,ll p){
	ll sum=0;
	while(x){
		if(x&1) sum=(sum+y)%p;
		y=(y*2)%p;
		x>>=1;
	}
	rt sum;
}
ll ksm(ll x,ll k,ll p){
	if(k==1) rt x;
	ll st=ksm(x,k/2,p);st=ksc(st,st,p);
	if(k&1) rt ksc(st,x,p);
	rt st;
}
I pd(ll x){
	if(x==2) rt 1;
	if(x<2||x%2==0) rt 0;
	ll d=x-1,y,k=0,s;
	while(d%2==0){k++,d/=2;}
	F(i,0,12){
		if(x==a[i]) rt 1;
		s=ksm(a[i],d,x);
		F(j,1,k){
			y=ksc(s,s,x);
			if(y==1&&s!=1&&s!=x-1) rt 0;
			s=y;
		}
		if(s!=1) rt 0;
	}
	rt 1;
}
I main(){
	P("prime");
	rd(n);
	while(n--){
		rd(x);
		if(pd(x)) printf("Prime\n");
		else printf("Not prime\n");
	}
	rt 0;
}

作者：zsjzliziyang
QQ:1634151125
转载及修改请注明
本文地址:https://blog.csdn.net/zsjzliziyang/article/details/99763123