两个有关素数的算法

Miller_Rabin算法（素数判定）

作用

单独判断一个大数是否素数。缺点他是一个不保证正确的算法，我们只能通过多次执行算法让这个错误的概率很小，不过幸运的是通常来看它的错误概率可以小到忽略不计。

补充:【模板】素数测试（Miller-Rabin测试）

算法的理论基础

Fermat定理：若n是奇素数，a是任意正整数(1≤ a≤ n−1)，则 a^(n-1) ≡ 1 mod n。

推演自Fermat定理, 如果n是一个奇素数，将n−1表示成2^s*r的形式，r是奇数，a与n是互素的任何随机整数，那么a^r ≡ 1 mod n或者对某个j (0 ≤ j≤ s−1, j∈Z) 等式a^(2jr) ≡ −1 mod n 成立。

代码实现

 1 bool Miller_Rabin(ll n)//Miller_Rabin算法，判断n是否为素数
 2 {
 3     if(n<2) return false;
 4     if(n==2) return true;
 5     if(!(n&1)) return false;
 6     ll r=n-1,s=0;
 7     while(!(r&1)){r=r>>1;s++;}
 8     for(ll i=0;i<NUM;i++)
 9     {
10         ll a=rand()%(n-1)+1;
11         if(check(a,n,r,s))
12             return false;
13     }
14     return true;
15 }

Pollard_rho算法

作用

求一个数的因子

原理

设n为待分解的大整数，用某种方法生成a和b，计算p=gcd(a-b,n),直到p不为1或a,b出现循环时为止，若p=n，则说明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。

可以在θ（sqrt(p)）的期望时间内找到n的一个小因子p。但对于因子很少，因子值很大的大整数n，该方法依然不是很有效。

为了减少反复的次数，算法做了一些改进。该算法用数对(x0,x0)开始，并且用x(i+1)=f(xi) ，迭代计算(x1,x2),(x2,x4),(x3,x6)…(xi,x2i)。在每一次迭代中，我们都应用上述函数式运算(从第二步)第一次计算数对中的第一个元素，第二次计算数对中的第二个元素。

代码实现

 1 ll Pollard_rho(ll n,ll c)//Pollard_rho算法，找出n的因子
 2 {
 3     ll i=1,j,k=2,d,p;
 4     ll x=rand()%n;
 5     ll y=x;
 6     while(true)
 7     {
 8         i++;
 9         x=(mul_mod(x,x,n)+c)%n;
10         if(y==x) return n;
11         if(y>x) p=y-x;
12         else p=x-y;
13         d=gcd(p,n);
14         if(d!=1&&d!=n) return d;
15         if(i==k)
16         {
17             y=x;
18             k+=k;
19         }
20     }
21 }
22 void find(ll n)//找出n的所有因子
23 {
24     if(Miller_Rabin(n))
25     {
26         f[t++]=n;//保存所有因子
27         return;
28     }
29     ll p=n;
30     while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);//由于p必定为合数，所以通过多次求解必定能求得答案
31     find(p);
32     find(n/p);
33 }

因数个数求解公式

对于任何一个自然数N，都可以分解质因子得到如下形式：
$N = p e 1 1 * p e 2 2 * p e 3 3 * \cdot \cdot \cdot * p e k k$
那么，N的因子的个数为：
$f (n) = (1 + e 1) * (1 + e 2) * \cdot \cdot \cdot * (1 + e k)$
如N=100，分解质因子变形为：100=22∗52，N的因子的个数为：f(N)=f(100)=(1+2)∗(1+2)=9。

即：1,2,4,5,10,20,25,50,100。

一些函数

头文件&宏定义

1 #include<bits/stdc++.h>
2 using namespace std;
3 #define ll long long
4 const ll NUM=10;//运算次数，Miller_Rabin算法为概率运算，误判率为2^(-NUM);
5 ll t,f[100];

快速乘法、快速幂函数

 1 ll mul_mod(ll a,ll b,ll n)//求a*b%n，由于a和b太大，需要用进位乘法
 2 {
 3     a=a%n;
 4     b=b%n;
 5     ll s=0;
 6     while(b)
 7     {
 8         if(b&1)
 9             s=(s+a)%n;
10         a=(a<<1)%n;
11         b=b>>1;
12     }
13     return s;
14 }
15 ll pow_mod(ll a,ll b,ll n)//求a^b%n
16 {
17     a=a%n;
18     ll s=1;
19     while(b)
20     {
21         if(b&1)
22             s=mul_mod(s,a,n);
23         a=mul_mod(a,a,n);
24         b=b>>1;
25     }
26     return s;
27 }
28 bool check(ll a,ll n,ll r,ll s)
29 {
30     ll ans=pow_mod(a,r,n);
31     ll p=ans;
32     for(ll i=1;i<=s;i++)
33     {
34         ans=mul_mod(ans,ans,n);
35         if(ans==1&&p!=1&&p!=n-1)
36             return true;
37         p=ans;
38     }
39     if(ans!=1) return true;
40     return false;
41 }
42 
43 ll gcd(ll a,ll b)
44 {
45     return b==0?a:gcd(b,a%b);
46 }

求因数个数时主函数

 1 int main()
 2 {
 3     srand(time(NULL));//随机数设定种子
 4     ll n;cin>>n;
 5     if(n==1){cout<<"1"<<endl;return 0;}
 6     t=0;
 7     find(n);
 8     sort(f,f+t);
 9     map<ll,int>q;
10     for(int i=0;i<t;i++)
11     {
12         q[f[i]]++;
13     }
14     map<ll,int>::iterator it;
15     ll ans=1;
16     for(it=q.begin();it!=q.end();it++)
17     {
18         int s=it->second;
19         ans*=1+s;
20     }
21     cout<<ans<<endl;
22     return 0;
23 }