大素数分解与素数判定

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. //噫  是不是狠有趣(。＾▽＾)  

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>//这个必须有
using namespace std;
int main()
{
    int a = 100;
    srand( time(NULL));
    while(a--)
        cout << rand() << endl;
    return 0;
}
//这个程序的作用是产生100个随机数
//如果你和我一样有颗童心去多试几次的话你会发现——每次产生的随机数都不一样
//噫  是不是狠有趣(。＾▽＾)

学了这么多是不是手痒了？别着急，点我有惊喜。

AC代码（C++【因为涉及到ctime，所以G++会RE的】）：

[cpp] view plain copy print ?

1. /* ************************************************* 
2.  * 
3.  * Miller_Rabin 算法进行素数测试  
4.  * 速度快，可以判断一个 < 2^63 的数是不是素数 
5.  * 
6.  **************************************************/  
7. #include <cstdio>  
8. #include <iostream>  
9. #include <cstring>  
10. #include <algorithm>  
11. #include <cmath>  
12. #include <ctime>  
13. using namespace std;  
14. const int S = 8;//随机算法判定次数，一般8~10次就够了  
15.   
16. //计算ret = (a*b)%c  a, b, c < 2^63  
17. long long mult_mod(long long a, long long b, long long c)  
18. {  
19.     a %= c;  
20.     b %= c;  
21.     long long ret = 0;  
22.     long long tem = a;  
23.     while(b)  
24.     {  
25.         if(b & 1)  
26.         {  
27.             ret += tem;  
28.             if(ret > c) ret -= c;//直接取模慢很多  
29.         }  
30.         tem <<= 1;  
31.         if(tem > c) tem -= c;  
32.         b >>= 1;  
33.     }  
34.     return ret;  
35. }  
36.   
37. //计算 ret = (a^n) % mod  
38. long long pow_mod(long long a, long long n, long long mod)  
39. {  
40.     long long ret = 1;  
41.     long long tem = a % mod;  
42.     while(n)  
43.     {  
44.         if(n & 1) ret = mult_mod(ret, tem, mod);  
45.         tem = mult_mod(tem, tem, mod);  
46.         n >>= 1;  
47.     }  
48.     return ret;  
49. }  
50.   
51. // 通过 a^(n-1)=1(mod n)来判断n是不是素数  
52. // n-1 = x * 2^t 中间使用二次判断  
53. // 是合数返回true，不一定是合数返回false  
54. bool check(long long a, long long n, long long x, long long t)  
55. {  
56.     long long ret = pow_mod(a, x, n);//a^x % n  
57.     long long last = ret;  
58.     for(int i = 1; i <= t; ++i)//进行t次(a^x % n)^2 % n  
59.     {  
60.         ret = mult_mod(ret, ret, n);  
61.         if(ret == 1 && last != 1 && last != n - 1) return true;//合数  
62.         last = ret;  
63.     }  
64.     if(ret != 1) return true;  
65.     return false;//不一定是合数  
66. }  
67.   
68. //**************************************************  
69. // Miller_Rabin算法  
70. // 是素数返回true,(可能是伪素数)  
71. // 不是素数返回false  
72. //**************************************************  
73. bool Miller_Rabin(long long n)  
74. {  
75.     if(n < 2) return false;  
76.     if(n == 2) return true;  
77.     if( (n&1) == 0) return false;//偶数  
78.     long long x = n - 1;  
79.     long long t = 0;  
80.     while( (x&1) == 0) //将n分解为x*2^t;  
81.     {  
82.         x >>= 1;  
83.         t++;  
84.     }  
85.   
86.     srand( time(NULL));  
87.   
88.     for(int i = 0; i < S; ++i)  
89.     {  
90.         long long a = rand()%(n-1) + 1;//产生随机数a(并控制其范围在1 ~ n-1之间)  
91.         if(check(a, n, x, t))//是合数  
92.             return false;  
93.     }  
94.     return true;  
95. }  
96.   
97. //**********************************************    
98. //    
99. // pollard_rho 算法进行质因素分解    
100. //    
101. //*********************************************   
102. int tol;//质因数的个数，编号为0~tol-1  
103. long long factor[100];//质因素分解结果(刚返回时是无序的)  
104.   
105. long long gcd(long long a, long long b)  
106. {  
107.     long long t;  
108.     while(b)  
109.     {  
110.         t = a;  
111.         a = b;  
112.         b = t % b;  
113.     }  
114.     if(a >= 0) return a;  
115.     return -a;  
116. }  
117.   
118. //找出一个因子  
119. long long pollard_rho(long long x, long long c)  
120. {  
121.     long long i = 1, k = 2;  
122.     srand( time(NULL));  
123.     long long x0 = rand()%(x-1) + 1;//产生随机数x0(并控制其范围在1 ~ x-1之间)  
124.     long long y = x0;  
125.     while(1)  
126.     {  
127.         i++;  
128.         x0 = (mult_mod(x0, x0, x) + c) % x;  
129.         long long d = gcd(y - x0, x);  
130.         if(d != 1 && d != x) return d;  
131.         if(y == x0) return x;  
132.         if(i == k)  
133.         {  
134.             y = x0;  
135.             k += k;  
136.         }  
137.     }  
138. }  
139.   
140. //对n进行素因子分解，存入factor。 k设置为107左右即可  
141. void findfac(long long n, int k)  
142. {  
143.     if(n == 1) return ;  
144.     if(Miller_Rabin(n))//是素数就把这个素因子存起来  
145.     {  
146.         factor[tol++] = n;  
147.         return ;  
148.     }  
149.     int c = k;  
150.     long long p = n;  
151.     while(p >= n)  
152.         p = pollard_rho(p, c--);//值变化，防止陷入死循环k  
153.     findfac(p, k);  
154.     findfac(n/p, k);  
155. }  
156.   
157. int main()  
158. {  
159.     int T;  
160.     long long n;  
161.     scanf("%d",&T);  
162.     while(T--)  
163.     {  
164.         scanf("%lld",&n);  
165.         if(Miller_Rabin(n)) cout << "Prime" << endl;  
166.         else  
167.         {  
168.             tol = 0;  
169.             findfac(n, 107);  
170.             long long ans = factor[0];  
171.             for(int i = 1; i < tol; ++i)  
172.                 ans = min(ans, factor[i]);  
173.             cout << ans << endl;  
174.         }  
175.     }  
176.     return 0;  
177. }  

/* *************************************************
 *
 * Miller_Rabin 算法进行素数测试 
 * 速度快，可以判断一个 < 2^63 的数是不是素数
 *
 **************************************************/
#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
using namespace std;
const int S = 8;//随机算法判定次数，一般8~10次就够了

//计算ret = (a*b)%c  a, b, c < 2^63
long long mult_mod(long long a, long long b, long long c)
{
    a %= c;
    b %= c;
    long long ret = 0;
    long long tem = a;
    while(b)
    {
        if(b & 1)
        {
            ret += tem;
            if(ret > c) ret -= c;//直接取模慢很多
        }
        tem <<= 1;
        if(tem > c) tem -= c;
        b >>= 1;
    }
    return ret;
}

//计算 ret = (a^n) % mod
long long pow_mod(long long a, long long n, long long mod)
{
    long long ret = 1;
    long long tem = a % mod;
    while(n)
    {
        if(n & 1) ret = mult_mod(ret, tem, mod);
        tem = mult_mod(tem, tem, mod);
        n >>= 1;
    }
    return ret;
}

// 通过 a^(n-1)=1(mod n)来判断n是不是素数
// n-1 = x * 2^t 中间使用二次判断
// 是合数返回true，不一定是合数返回false
bool check(long long a, long long n, long long x, long long t)
{
    long long ret = pow_mod(a, x, n);//a^x % n
    long long last = ret;
    for(int i = 1; i <= t; ++i)//进行t次(a^x % n)^2 % n
    {
        ret = mult_mod(ret, ret, n);
        if(ret == 1 && last != 1 && last != n - 1) return true;//合数
        last = ret;
    }
    if(ret != 1) return true;
    return false;//不一定是合数
}

//**************************************************
// Miller_Rabin算法
// 是素数返回true,(可能是伪素数)
// 不是素数返回false
//**************************************************
bool Miller_Rabin(long long n)
{
    if(n < 2) return false;
    if(n == 2) return true;
    if( (n&1) == 0) return false;//偶数
    long long x = n - 1;
    long long t = 0;
    while( (x&1) == 0) //将n分解为x*2^t;
    {
        x >>= 1;
        t++;
    }

    srand( time(NULL));

    for(int i = 0; i < S; ++i)
    {
        long long a = rand()%(n-1) + 1;//产生随机数a(并控制其范围在1 ~ n-1之间)
        if(check(a, n, x, t))//是合数
            return false;
    }
    return true;
}

//**********************************************  
//  
// pollard_rho 算法进行质因素分解  
//  
//********************************************* 
int tol;//质因数的个数，编号为0~tol-1
long long factor[100];//质因素分解结果(刚返回时是无序的)

long long gcd(long long a, long long b)
{
    long long t;
    while(b)
    {
        t = a;
        a = b;
        b = t % b;
    }
    if(a >= 0) return a;
    return -a;
}

//找出一个因子
long long pollard_rho(long long x, long long c)
{
    long long i = 1, k = 2;
    srand( time(NULL));
    long long x0 = rand()%(x-1) + 1;//产生随机数x0(并控制其范围在1 ~ x-1之间)
    long long y = x0;
    while(1)
    {
        i++;
        x0 = (mult_mod(x0, x0, x) + c) % x;
        long long d = gcd(y - x0, x);
        if(d != 1 && d != x) return d;
        if(y == x0) return x;
        if(i == k)
        {
            y = x0;
            k += k;
        }
    }
}

//对n进行素因子分解，存入factor。 k设置为107左右即可
void findfac(long long n, int k)
{
    if(n == 1) return ;
    if(Miller_Rabin(n))//是素数就把这个素因子存起来
    {
        factor[tol++] = n;
        return ;
    }
    int c = k;
    long long p = n;
    while(p >= n)
        p = pollard_rho(p, c--);//值变化，防止陷入死循环k
    findfac(p, k);
    findfac(n/p, k);
}

int main()
{
    int T;
    long long n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld",&n);
        if(Miller_Rabin(n)) cout << "Prime" << endl;
        else
        {
            tol = 0;
            findfac(n, 107);
            long long ans = factor[0];
            for(int i = 1; i < tol; ++i)
                ans = min(ans, factor[i]);
            cout << ans << endl;
        }
    }
    return 0;
}

经过了“几天”的学习，终于能明白M-r,p-r算法是怎么实现的了（吐槽一下那三位对我关爱有加的兄弟，给我留了个这么有用的知识点让我讲），然后对着板子敲了几遍熟悉了一下。可能是还没碰到这种题目吧，内心里总觉得。。？总之也算是没浪费这些时间，至少我可以对着板子来对这知识点进行自在应用了（比如用朴素算法一不小心就会超时的一道题目。