两个有关素数的算法
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=pe11∗pe22∗pe33∗···∗pekkN=p1e1∗p2e2∗p3e3∗···∗pkek那么,N的因子的个数为:
f(n)=(1+e1)∗(1+e2)∗···∗(1+ek)f(n)=(1+e1)∗(1+e2)∗···∗(1+ek)如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 }
文章搬运自我的个人博客http://duny31030.top
原博客为静态博客,因备份丢失无法继续更新,所以又搬运回博客园,可能部分文章阅读体验不好,可以到我的静态博客搜索相同标题查看