给你一个数N(2 <= N < 2^54) ,若N是素数 , 输出Prime, 否则输出最小的素因子。
由于N很大,所以只能先用Miller_Rabin算法进行素数判断,然后用Pollard_rho分解因子。
Miller-Rabin算法本质上是一种概率算法,存在误判的可能性,但是出错的概率非常小。出错的概率到底是多少,存在严格的理论推导。
费马小定理
- 如果p是质数且(a,p)=1,则有a^p−1≡1(mod p)。
当然反过来不一定成立。即当(a^p−1)%p=1时,p未必是质数。但是这个概率比较小。所以利用费尔马小定理来检测素数,不能保证时刻都对,只能保证出错的概率比较小。
给定正整数n,问n是否为质数(显然只需判断正奇数),最基本的做法就是计算2^(n−1)%n是否为1。如果不是1,n肯定为合数;否则,n可能为质数。
二次探测定理
- 如果p是一个素数且x^2≡1(mod p) , 则 x=1或者x=−1,注意到在模p的意义下,x=−1等价于x=p−1。
- 证明:方程x^2≡1(mod p)可以改写为(x + 1)(x - 1) % p = 0.所以x为1或-1.但在mod p的情况下,x=-1 等价于 x=p-1
Miller-Rabin算法
利用上面两个定理,就可以构造出Miller-Rabin算法。考虑到n肯定是奇数,则n一定可以表示为n−1=2^s∗d,其中s≥1且d是奇数。则
a^(n−1)=a^(2^s∗d)=(((a^d)^2)^...^)2
也就是说,a^(n−1)相当于a^d平方若干次。例如当n=7时,a^(n−1)就是a^6,就是a^3的平方。当n=13时,a^(n−1)就是a^12,就是a^3的平方的平方。
以n=13的情况进行说明(所有运算都是在模n的意义下,以下的文字说明省略了这一点),任取一个a , 1< a <13,计算a^3,再将其平方一次得到a^6,注意到a^3是a^6的平方根,根据平方根定理的推论,如果a^6=1且a^3≠±1,则n肯定是合数。将a^6平方一次得到a^12,同样,如果a^12=1且a^6≠±1,则n肯定是合数。最后,根据费尔马小定理,如果a^12≠1,则n肯定是合数。否则,n有极大概率为质数。
为了增加得到正确判断的概率,可以将a重复取不同的值,对每一个a验证一次a^d到a^(n−1)的过程。不过,考虑到ACM的特殊性,测试数据应该不会选择伪素数特别是强伪素数。所以很多题目的AC程序实际上只来一次即可。
代码
typedef long long ll;
//利用二进制计算a*b%mod
ll multiMod(ll a,ll b,ll mod){
a %= mod , b %= mod;
ll ret = 0;
while(b){
if (b & 1) {
ret += a;
if(ret >= mod) ret -= mod;
}
a <<= 1;
if(a >= mod) a -= mod;
b >>= 1;
}
return ret;
}
//计算a^b%mod
ll powerMod(ll a,ll b,ll mod){
ll ret = 1;
a %= mod;
while(b){
if (b & 1) ret = multiMod(ret , a , mod);
b >>= 1;
a = multiMod(a , a , mod);
}
return ret;
}
//Miller-Rabin测试,测试n是否为素数
bool Miller_Rabin(ll n,int repeat){
if (2 == n || 3 == n) return true;
if (!(n & 1)) return false;
//将n分解为2^s*d
ll d = n - 1;
int s = 0;
while((d&1)== 0) ++s, d>>=1LL;
//srand((unsigned)time(0)); G++不能使用srand(time(NULL)) , 会RE
for(int i = 0; i < repeat; ++i){//重复repeat次
ll a = rand() % (n - 3) + 2;//取一个随机数,[2,n-1)
ll x = powerMod(a , d , n);
ll y = 0;
for(int j = 0; j < s; ++j){
y = multiMod(x , x , n);
if ( 1 == y && 1 != x && n-1 != x ) return false;
x = y;
}
if ( 1 != y ) return false;
}
return true;
}
参考:http://blog.csdn.net/u012061345/article/details/48241581#
pollard_rho 算法
- 通过某种方法得到两个整数a和b,而待分解的大整数为n,计算p=gcd(a-b,n),直到p不为1,或者a,b出现循环为止。然后再判断p是否为n,如果p=n成立,那么返回n是一个质数,否则返回p是n的一个因子,那么我们又可以递归的计算Pollard(p)和Pollard(n/p),这样,我们就可以求出n的所有质因子。
- 具体操作中,我们通常使用函数x2=x1*x1+c来计算逐步迭代计算a和b的值,实践中,通常取c为1,即b=a*a+1,在下一次计算中,将b的值赋给a,再次使用上式来计算新的b的值,当a,b出现循环时,即可退出进行判断。
- 在实际计算中,a和b的值最终肯定一出现一个循环,而将这些值用光滑的曲线连接起来的话,可以近似的看成是一个ρ型的。
对于Pollard rho,它可以在O(sqrt(p))的时间复杂度内找到n的一个小因子p,可见效率还是可以的,但是对于一个因子很少、因子值很大的大整数n来说,Pollard rho算法的效率仍然不是很好,那么,我们还得寻找更加的方法了。
代码
ll factor[100]; //素因子的结果,返回时是无序的
int tol; //素因子的个数
ll gcd(ll a , ll b)
{
ll t;
while(b) {t = a; a = b; b = t%b;}
if(a >= 0) return a;
return -a;
}
ll pollard_rho(ll n , ll c)
{
ll i = 1 , k = 2;
//srand(time(NULL));
ll a = rand() % (n - 1) + 1 , b = a;
while(1)
{
i++;
a = (multiMod(a , a , n) + c) % n;
ll p = gcd(b - a , n);
if(p != 1 && p != n) return p;
if(b == a) return n;
if(i == k) {b = a; k += k;}
}
}
void findfac(ll n , int k)//k为自己设定的值,一般取1
{
if(1 == n) return ;
if(Miller_Rabin(n , 8))
{
factor[tol++] = n;
return ;
}
ll p = n;
int c = k;
while(p >= n) p = pollard_rho(p , c--);
findfac(p , k);
findfac(n/p , k);
}