Miller-Rabin算法 codevs 1702 素数判定 2

转载自:http://www.dxmtb.com/blog/miller-rabbin/

普通的素数测试我们有O(√ n)的试除算法。事实上,我们有O(slog³n)的算法。

定理一:假如p是质数,且(a,p)=1,那么a^(p-1)≡1(mod p)。即假如p是质数,且a,p互质,那么a的(p-1)次方除以p的余数恒等于1。(费马小定理)

该定理的逆命题是不一定成立的,但是令人可喜的是大多数情况是成立的。

于是我们就得到了一个定理的直接应用,对于待验证的数p,我们不断取a∈[1,p-1]且a∈Z,验证a^(p-1) mod p是否等于1,不是则p果断不是素数,共取s次。其中a^(p-1) mod p可以通过把p-1写成二进制,由(a*b)mod c=(a mod c)*b mod c,可以在t=log(p-1)的时间内计算出解,如考虑整数相乘的复杂度,则一次计算的总复杂度为log³(p-1)。这个方法叫快速幂取模。

为了提高算法的准确性,我们又有一个可以利用的定理。
定理二:对于0<x<p,x^2 mod p =1 => x=1或p-1。

我们令p-1=(2^t)*u,即p-1为u二进制表示后面跟t个0。我们先计算出x[0]=a^u mod p ,再平方t次并在每一次模p,每一次的结果记为x[i],最后也可以计算出a^(p-1) mod p。若发现x[i]=1而x[i-1]不等于1也不等于p-1,则发现p果断不是素数。

可以证明,使用以上两个定理以后,检验s次出错的概率至多为2^(-s),所以这个算法是很可靠的。

需要注意的是,为了防止溢出(特别大的数据),a*b mod c 也应用类似快速幂取模的方法计算。当然,数据不是很大就可以免了

 

codevs 1702 素数判定 2

 时间限制: 1 s
 空间限制: 128000 KB
 题目等级 : 钻石 Diamond
 
题目描述  Description

一个数,他是素数么?

设他为P满足(P<=2^63-1)

 

输入描述  Input Description

P

输出描述  Output Description

Yes|No

样例输入  Sample Input

2

 

样例输出  Sample Output

Yes

 

数据范围及提示  Data Size & Hint

算法导论——数论那一节
注意Carmichael Number

分类标签 Tags 点此展开 
 素数判定 数论
 
解析:因为数据范围是2^63-1很大,一般方法时间复杂度很高,所以用随机化算法--Miller Rabin算法
 1 /*数据确实没有错,我的错误点是,在两个long long相乘的时候,直接进行了乘法运算,很有可能溢出,所以错了后面的点,long long数据要再写一个相乘%的子函数*/
 2 #include<iostream>
 3 using namespace std;
 4 #include<cstdio>
 5 #define S 10
 6 #include<cstdlib>
 7 #include<ctime>
 8 #define ll long long
 9 ll cas, maxz;
10 ll read()
11 {
12     ll ans=0;char c;
13     c=getchar();
14     while(c<'0'||c>'9') c=getchar();
15     while(c>='0'&&c<='9') 
16     {
17         ans=ans*10+c-'0';
18         c=getchar();
19     }
20     return ans;
21 }
22 ll quick_mul_mod(ll a,ll b,ll c)//a*b%c
23 {
24     ll ret=0;
25     a%=c;b%=c;
26     while(b)
27     {
28         if(b&1)
29         {
30             ret+=a;/*相当于加了b次a*/
31             ret%=c;
32             b--;
33         }
34         a<<=1;
35         a%=c;
36         b>>=1;
37     }
38     return ret;
39 }
40 ll quick_mod(ll a,ll b,ll c)//ji suan a^b%c
41 {
42     ll ans=1;
43     a%=c;
44     while(b)
45     {
46         if(b&1)
47         {
48             b--;
49             ans=quick_mul_mod(ans,a,c);/*注意只要是long long的乘法,很有可能溢出的*/
50         }
51         b>>=1;
52         a=quick_mul_mod(a,a,c);
53     }
54     return ans;
55 }
56 bool Miller_rabin(ll n)
57 {
58     if(n==2) return true;
59     if(n<=1||!(n&1)) return false;
60     ll u=n-1,t=0;
61     while(!(u&1))
62     {
63         u>>=1;
64         t++;
65     }
66     for(int i=0;i<S;++i)
67     {
68         ll x=rand()%(n-1)+1;
69         x=quick_mod(x,u,n);
70         for(int i=1;i<=t;++i)
71         {
72             ll y=quick_mul_mod(x,x,n);
73             if(y==1&&x!=1&&x!=n-1)
74               return false;
75             x=y;
76         }
77         if(x!=1) return false;
78     }
79     return true;
80 }
81 
82 int main()
83 {
84         ll n=read();
85         if(Miller_rabin(n))
86           printf("Yes\n");
87         else printf("No\n");
88     return 0;
89  }

 

转载于:https://www.cnblogs.com/c1299401227/p/5514696.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值