[模板] Miller_Rabin和Pollard_Rho

Miller_Rabin

快速($O(slogn)$,s为尝试次数)地测试一个数是否是质数

 

首先有费马小定理$a^{p-1}=1\ (mod\ p)$当p为质数时成立,所以可以随机选择a来以这个式子作为一定的判断依据,但并不是所有合数都不满足这个式子,甚至存在合数对所有的a都不满足这个式子

然后有二次探测定理:若p是质数,$a^2=1\ (mod\ p)$成立当且仅当$a=1\ (mod\ p)$或$a=p-1\ (mod\ p)$

证明:移项可得$(a-1)(a+1)=0\ (mod\ p)$,于是有$p|(a-1)(a+1)$,若p为质数,则只能有$p|(a-1)$或者$p(a+1)$

所以,如果有$a^2=1\ (mod\ p)$但a不等于1或者p-1的话,证明p一定是合数

然而当我们随便选一个a,很难保证它的平方就等于1...但我们还有费马小定理!

所以我们可以考虑把p-1表示成$s*2^t$的形式,然后从$a^s$一路平方上去,并进行判断

这个a可以考虑选择前10个质数,在1e18的范围内应该都是对的

 

Pollard_Rho

快速(期望$O(n^{1/4})$)地找到n的某个非平凡因子

 

首先,用Miller_Rabin来判断n是否是质数(这也是P_R和M_R唯一的关联)

然后我们考虑随机一个数x,然后计算$gcd(x,n)$,如果结果不是1或者n,就找到了n的一个非平凡因子

然而纯随机肯定是不行的..

考虑一种生成随机数的方式:$x_i=x_{i-1}^2+c\ (mod\ n)$,$c$和$x_0$为随机的常数,这样生成出的数一定会成一个"ρ"形,而且根据生日悖论(期望随机$\sqrt{n}$次就能得到在n的范围内相等的两个数),环长和前面一段的长度都是期望$\sqrt{n}$的

更进一步的是,如果考虑a是n的一个最小质因数,那么$x_i$在模a意义下,会形成一个期望$\sqrt{a}$的ρ

这样的话,我随机在这个模n的ρ形上取两个数做gcd,期望只要$\sqrt{a}$也就是$n^{1/4}$次,就能找到一个n的因数

于是可以使用a走一步,b走两步的方法判圈,同时计算$gcd(a-b,n)$

但注意到一次gcd需要log的时间,所以我们可以采用每127(我也不知道这个玄学的数是哪来的)个计算一次gcd的方法来优化,就是说,每次连续算127个$a-b$并乘起来,然后算一下gcd,如果发现gcd大于1了,再重新算一下这127个看看是哪个大于1了

 

例题

luogu4718 (求n的最大质因子)
 1 #include<bits/stdc++.h>
 2 #include<tr1/unordered_map>
 3 #define CLR(a,x) memset(a,x,sizeof(a))
 4 #define MP make_pair
 5 #define fi first
 6 #define se second
 7 using namespace std;
 8 typedef long long ll;
 9 typedef unsigned long long ull;
10 typedef long double ld;
11 typedef pair<int,int> pa;
12 const int maxn=233;
13 
14 inline ll rd(){
15     ll x=0;char c=getchar();int neg=1;
16     while(c<'0'||c>'9'){if(c=='-') neg=-1;c=getchar();}
17     while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar();
18     return x*neg;
19 }
20 
21 int pri[10]={2,3,5,7,11,13,17,19,23,29};
22 
23 inline ll mul(ll a,ll b,ll p){
24     return ((a*b-(ll)((ld)a*b/p)*p)%p+p)%p;
25 }
26 
27 inline ll fpow(ll a,ll b,ll p){
28     ll r=1;
29     while(b){
30         if(b&1) r=mul(r,a,p);
31         a=mul(a,a,p),b>>=1;
32     }return r;
33 }
34 
35 inline bool judge(ll p){
36     if(p==2) return 1;
37     else if(p==1||p%2==0) return 0;
38     ll s=0,t=p-1;
39     while(!(t&1)) s++,t>>=1;
40     for(int i=0;i<10&&pri[i]<p;i++){
41         ll x=fpow(pri[i],t,p);
42         for(int j=1;j<=s;j++){
43             ll r=mul(x,x,p);
44             if(r==1&&x!=1&&x!=p-1) return 0;
45             x=r;
46         }
47         if(x!=1) return 0;
48     }return 1;
49 }
50 
51 inline ll gcd(ll a,ll b){return a?gcd(b%a,a):b;}
52 inline ll rand(ll r){
53     return (((ll)rand()<<30)|rand())%(r-1)+1;
54 }
55 
56 inline ll calc(ll x){
57     if(judge(x)) return x;
58     while(1){
59         ll a=rand(x-1),c=rand(x-1),b=(mul(a,a,x)+c)%x;
60         // printf("~%lld %lld %lld\n",a,b,c);
61         bool bl=0;
62         while(a!=b){
63             ll s=1;ll aa=a,bb=b;
64             if(!bl){
65                 for(int i=1;i<=127&&aa!=bb;i++){
66                     s=mul(s,aa-bb,x);
67                     aa=(mul(aa,aa,x)+c)%x;
68                     bb=(mul(bb,bb,x)+c)%x,bb=(mul(bb,bb,x)+c)%x;
69                 }
70             }
71             ll g=gcd(s+x,x);
72             if(g>1) bl=1;
73             else a=aa,b=bb;
74             if(bl){
75                 g=gcd(a-b+x,x);
76                 if(g>1) return max(calc(g),calc(x/g));
77                 a=(mul(a,a,x)+c)%x;
78                 b=(mul(b,b,x)+c)%x,b=(mul(b,b,x)+c)%x;
79             }
80             
81             
82         }
83     }
84 }
85 
86 int main(){
87     srand(1919810);
88     for(int T=rd();T;T--){
89         ll x=rd();
90         ll re=calc(x);
91         if(re==x) puts("Prime");
92         else printf("%lld\n",re);
93     }
94     return 0;
95 }

 

posted @ 2019-06-24 13:33  Ressed  阅读(109)  评论(0编辑  收藏