素數判定與大數分解【Miller-rabin算法】【pollard-rho算法】

對應練習題：SDNUOJ1 1286

點選打開連結

1.Miller-rabin算法：

Miller-rabin算法是一個用來快速判斷一個正整數是否為素數的算法。

根據費馬小定理，如果p是素數，則a^(p-1)≡1(mod p)對所有的a∈[1,n-1]成立。是以如果在[1,n-1]中随機取出一個a，發現不滿足費馬小定理，則證明n必為合數。

【但是每次嘗試過程中還做了一個優化操作，以提高用少量的a檢測出p不是素數的機率。這個優化叫做二次探測。它是根據這個定理：如果p是一個素數，那麼對于x(0<x<p)，若x^2%p=1，則x=1或p-1。】

為了計算a^(n-1)mod n，我們把n-1分解為x* 2^t的形式，其中t>=1且x是奇數；是以，a^(n-1)≡(a^x)^(2^t)(mod n),是以可以通過先計算a^x mod n,然後對結果連續平方t次來計算a^(n-1) mod n。一旦發現某次平方後mod n等于1了，那麼說明符合了二次探測定理的逆否命題使用條件，立即檢查x是否等于1或n-1，如果不等于1也不等于n-1則可直接判定p為合數。

2.pollard-rho算法：

這是一個用來快速對整數進行質因數分解的算法，需要與Miller-rabin共同使用。

算法原理：

1.通過某種方法得到兩個整數a和b，而待分解的大整數為n。

2.計算p=gcd(a-b,n)，直到p不為1(就是a-b與n不是互質)，或者a，b出現循環為止。

3.然後再判斷p=n？

4.如果p=n，那麼傳回n是一個質數。

5.否則傳回p是n的一個因子，那麼我們又可以遞歸的計算Pollard(p)和Pollard(n/p)，這樣，我們就可以求出n的所有質因子。

算法步驟：選取一個小的随機數x1，疊代生成x[i] = x[i-1]^2+c，一般取c=1，若序列出現循環則退出，計算p=gcd(x[i-1]-x[i],n)，若p=1則傳回上一步繼續疊代，否則跳出疊代過程。若p=n，則n為素數，否則p為n的一個約數，并遞歸分解p和n/p。

【小知識】：随機數生成

C++中函數srand（），可以指定不同的數（無符号整數變元）為種子。但是如果種子相同，僞随機數列也相同。

比較理想的是用變化的數，比如時間來作為随機數生成器的種子。 time的值每時每刻都不同，即種子不同，是以，産生的随機數也不同。 

用法什麼的想深入了解自己去搜吧，這裡隻要明白下面的程式中随機數是這樣産生的就行了。然後，在這裡再舉個小栗子以加深一下對它的了解：

[cpp]  view plain  copy  print ?

1. #include <cstdio>  
2. #include <iostream>  
3. #include <cstring>  
4. #include <algorithm>  
5. #include <cmath>  
6. #include <ctime>//這個必須有  
7. using namespace std;  
8. int main()  
9. {  
10.     int a = 100;  
11.     srand( time(NULL));  
12.     while(a--)  
13.         cout << rand() << endl;  
14.     return 0;  
15. }  
16. //這個程式的作用是産生100個随機數  
17. //如果你和我一樣有顆童心去多試幾次的話你會發現——每次産生的随機數都不一樣  
18. //噫  是不是狠有趣(。＾▽＾)  

學了這麼多是不是手癢了？别着急，點我有驚喜。

AC代碼（C++【因為涉及到ctime，是以G++會RE的】）：

[cpp]  view plain  copy  print ?

1. #include <cstdio>  
2. #include <iostream>  
3. #include <cstring>  
4. #include <algorithm>  
5. #include <cmath>  
6. #include <ctime>  
7. using namespace std;  
8. const int S = 8;//随機算法判定次數，一般8~10次就夠了  
9. //計算ret = (a*b)%c  a, b, c < 2^63  
10. long long mult_mod(long long a, long long b, long long c)  
11. {  
12.     a %= c;  
13.     b %= c;  
14.     long long ret = 0;  
15.     long long tem = a;  
16.     while(b)  
17.     {  
18.         if(b & 1)  
19.         {  
20.             ret += tem;  
21.             if(ret > c) ret -= c;//直接取模慢很多  
22.         }  
23.         tem <<= 1;  
24.         if(tem > c) tem -= c;  
25.         b >>= 1;  
26.     }  
27.     return ret;  
28. }  
29. //計算 ret = (a^n) % mod  
30. long long pow_mod(long long a, long long n, long long mod)  
31. {  
32.     long long ret = 1;  
33.     long long tem = a % mod;  
34.     while(n)  
35.     {  
36.         if(n & 1) ret = mult_mod(ret, tem, mod);  
37.         tem = mult_mod(tem, tem, mod);  
38.         n >>= 1;  
39.     }  
40.     return ret;  
41. }  
42. // 通過 a^(n-1)=1(mod n)來判斷n是不是素數  
43. // n-1 = x * 2^t 中間使用二次判斷  
44. // 是合數傳回true，不一定是合數傳回false  
45. bool check(long long a, long long n, long long x, long long t)  
46. {  
47.     long long ret = pow_mod(a, x, n);//a^x % n  
48.     long long last = ret;  
49.     for(int i = 1; i <= t; ++i)//進行t次(a^x % n)^2 % n  
50.     {  
51.         ret = mult_mod(ret, ret, n);  
52.         if(ret == 1 && last != 1 && last != n - 1) return true;//合數  
53.         last = ret;  
54.     }  
55.     if(ret != 1) return true;  
56.     return false;//不一定是合數  
57. }  
58. //**************************************************  
59. // Miller_Rabin算法  
60. // 是素數傳回true,(可能是僞素數)  
61. // 不是素數傳回false  
62. //**************************************************  
63. bool Miller_Rabin(long long n)  
64. {  
65.     if(n < 2) return false;  
66.     if(n == 2) return true;  
67.     if( (n&1) == 0) return false;//偶數  
68.     long long x = n - 1;  
69.     long long t = 0;  
70.     while( (x&1) == 0) //将n分解為x*2^t;  
71.     {  
72.         x >>= 1;  
73.         t++;  
74.     }  
75.     srand( time(NULL));  
76.     for(int i = 0; i < S; ++i)  
77.     {  
78.         long long a = rand()%(n-1) + 1;//産生随機數a(并控制其範圍在1 ~ n-1之間)  
79.         if(check(a, n, x, t))//是合數  
80.             return false;  
81.     }  
82.     return true;  
83. }  
84. //**********************************************    
85. //    
86. // pollard_rho 算法進行質因素分解    
87. //    
88. //*********************************************   
89. int tol;//質因數的個數，編号為0~tol-1  
90. long long factor[100];//質因素分解結果(剛傳回時是無序的)  
91. long long gcd(long long a, long long b)  
92. {  
93.     long long t;  
94.     while(b)  
95.     {  
96.         t = a;  
97.         a = b;  
98.         b = t % b;  
99.     }  
100.     if(a >= 0) return a;  
101.     return -a;  
102. }  
103. //找出一個因子  
104. long long pollard_rho(long long x, long long c)  
105. {  
106.     long long i = 1, k = 2;  
107.     srand( time(NULL));  
108.     long long x0 = rand()%(x-1) + 1;//産生随機數x0(并控制其範圍在1 ~ x-1之間)  
109.     long long y = x0;  
110.     while(1)  
111.     {  
112.         i++;  
113.         x0 = (mult_mod(x0, x0, x) + c) % x;  
114.         long long d = gcd(y - x0, x);  
115.         if(d != 1 && d != x) return d;  
116.         if(y == x0) return x;  
117.         if(i == k)  
118.         {  
119.             y = x0;  
120.             k += k;  
121.         }  
122.     }  
123. }  
124. //對n進行素因子分解，存入factor。 k設定為107左右即可  
125. void findfac(long long n, int k)  
126. {  
127.     if(n == 1) return ;  
128.     if(Miller_Rabin(n))//是素數就把這個素因子存起來  
129.     {  
130.         factor[tol++] = n;  
131.         return ;  
132.     }  
133.     int c = k;  
134.     long long p = n;  
135.     while(p >= n)  
136.         p = pollard_rho(p, c--);//值變化，防止陷入死循環k  
137.     findfac(p, k);  
138.     findfac(n/p, k);  
139. }  
140. int main()  
141. {  
142.     int T;  
143.     long long n;  
144.     scanf("%d",&T);  
145.     while(T--)  
146.     {  
147.         scanf("%lld",&n);  
148.         if(Miller_Rabin(n)) cout << "Prime" << endl;  
149.         else  
150.         {  
151.             tol = 0;  
152.             findfac(n, 107);  
153.             long long ans = factor[0];  
154.             for(int i = 1; i < tol; ++i)  
155.                 ans = min(ans, factor[i]);  
156.             cout << ans << endl;  
157.         }  
158.     }  
159.     return 0;  
160. }