Pollard_rho大数因数分解

从今天终于把这个算法整理出来了,比较弱理解起来有点慢(目前也不怎么懂),但是能愉快的运行了hhhh

题目链接https://nanti.jisuanke.com/t/28395 

Pollard_rho算法的大致流程是 先判断当前数是否是素数(Miller_rabin)了,如果是则直接返回。如果不是素数的话,试图找到当前数的一个因子(可以不是质因子)。然后递归对该因子和约去这个因子的另一个因子进行分解。

 

对于一个大整数n,采用一种随机化算法,我们取任意一个数x使得x是n的质因数的几率很小,但如果取两个数x1以及x2使得它们的差是n的因数的几率就提高了,如果取x1以及x2使得gcd(abs(x1−x2),n)>1的概率就更高了。这就是Pollard-Rho算法的主要思想。我们假设要找的因子为p,他是随机取一个x1,由x1构造x2,使得{p可以整除x1-x2 && x1-x2不能整除n}则p=gcd(x1-x2,n),结果可能是1也可能不是1。如果不是1就找寻成功了一个因子,返回因子;如果是1就寻找失败,那么我们就要不断调整x2,具体的办法通常是x2=x2*x2+c(c是自己定的)直到出现x2出现了循环==x1了表示x1选取失败重新选取x1重复上述过程。(似乎还存在一个每次找寻范围*2的优化,但是不太懂。。。)

那么我们怎样不断取得x1和x2呢? 
x1在区间[1,n]中随机(rand)出来,而x2则由x2=(x1*x1%n+c)%n推算出来,其中c为任意给定值,事实证明,这样就是比较优的。因为x1和x2再调整时最终一定会出现循环,形成一个类似希腊字母rho的形状,故因此得名。

 

 

另外通过factorization函数来分解素数,如果找到了一个素数因子则加入到因子map中,否则如果用Pollard找到一个因子则递归去找素数因子,具体解释可以看这个

代码如下:

#include <iostream> 
#include <algorithm>  
#include <cstring>  
#include <map> 
using namespace std;

typedef long long ll;
const int maxn = 100001;
const int times = 10;
int prime[] = { 2, 3, 5, 7, 13, 17, 11, 61, 97, 24251, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 67, 71, 73, 79, 83, 89 };		//得到随机检验算子 a
map<ll, int>ump;

ll  Random(ll n) {
	return (double(rand()) / RAND_MAX * n + 0.5);
}

ll q_mul(ll a, ll b, ll mod) {
	ll ans = 0;
	while (b) {
		if (b & 1)	ans = (ans + a) % mod;
		a = (a << 1) % mod;
		b >>= 1;
	}
	return ans;
	//	return (a*b)%mod;//超时就换成这个 
}

ll q_pow(ll a, ll b, ll mod) {
	ll ans = 1;
	while (b) {
		if (b & 1)	ans = q_mul(ans, a, mod);
		a = q_mul(a, a, mod);
		b >>= 1;
	}
	return ans;
}

bool miller_rabin(ll n) {				 //检验n是否是素数
	if (n<2) return false;
	if (n == 2) return true;
	if (!(n & 1)) return false;				//如果是2则是素数,如果<2或者是>2的偶数则不是素数

	ll k, j = 0, x, tem = n - 1;
	while (!(tem & 1)) tem >>= 1, j++;
	for (int i = 0; i<times; i++) {				//做times次随机检验
		if (prime[i] >= n) return true;
		x = q_pow(prime[i], tem, n);				//得到a^r mod n
		if (x == 1) continue;			//余数为1则为素数
		for (k = 0; k<j; k++) {
			if (x == n - 1) break;			//否则试验条件2看是否有满足的 j
			x = q_mul(x, x, n);
		}
		if (k == j) return false;
	}
	return true;
}

ll gcd(ll a, ll b) {
	if (a < 0) return gcd(-a, b);//加上这个能快些 
	return b ? gcd(b, a%b) : a;
}

ll pollard_rho(ll n, ll c) {//找到n的一个因子
	ll x, y, d, i = 1, k = 2;
	y = x = Random(n - 1) + 1;
	while (1) {
		i++;
		x = (q_mul(x, x, n) + c) % n;
		d = gcd(y - x, n);
		if (1<d&&d<n) return d;
		if (y == x) return n;//找到循环,选取失败,重新来
		if (i == k) {//似乎是一个优化,但是不是很清楚
			y = x;
			k <<= 1;
		}
	}
}

void factorization(ll n) {//分解 
	if (n == 1) return;		//此题特殊判断 
	if (miller_rabin(n)) {
		ump[n]++;
		return;
	}
	ll p = n;
	while (p >= n) p = pollard_rho(p, Random(n) + 1);
	factorization(p);
	factorization(n / p);
}

int main() {
	ll x, ret;
	while (cin >> x) {
		ret = 1;
		ump.clear();
		factorization(x);
		for (map<ll, int>::iterator it = ump.begin(); it != ump.end(); it++) {
			ret *= it->second + 1;
		}
		cout << ret << endl;
	}
	return 0;
}

参考:https://blog.csdn.net/maxichu/article/details/45459533

           https://blog.csdn.net/dllpxfire/article/details/81037933

  • 2
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值