1. 程式人生 > >Pollard-Rho演算法--大數分解

Pollard-Rho演算法--大數分解

其實剛開始知道這個演算法時,我以為需要字串操作什麼的,畢竟是大數嘛,可這傢伙只用了個long long,無語了....

long long int,能叫大數嗎,連2^100次方都處理不了。

說是這樣說,這個演算法已經不錯了,複雜度為O(n^1/4), 貌似目前學界沒有找到特別好的演算法,據說有什麼艾-魯法,威廉斯夫法,可以分解一個千位素數(需要一週時間)(上面兩種演算法百度不到QAQ),量子的那個shore法應該很快。

對於一個大整數n,我們取任意一個數xx使得xxnn的質因數的機率很小,但如果取兩個數x1x1以及x2

使得它們的差是n的因數的機率就提高了,如果x1以及x2使得gcd(abs(x1-x2), n) > 1的概率就更高了,這就是Pollard-Rho演算法

的思想。(概率的增加是因為組合數增加了)

又因為是質因數分解,所以我們要用到Miller_Rabbin演算法來判斷是否為素數,而隨機地取x1, x2,我們則是在[1~n]中隨機取x1,x2由x[i] = (x[i-1]*x[i-1]%n+c)%n,推算出來。

一個普遍模板如下:

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<iostream>
#include<queue>
using namespace std;

typedef long long ll;
const int maxn = 105;
ll x[maxn], ans;
queue<ll> aria;
ll min(ll a, ll b)
{
	if(a < b)	return a;
	else	return b;
}
ll multi(ll a, ll b, ll p)	//快速乘? 
{
	ll ans = 0;
	while(b){
		if(b & 1LL)	ans = (ans+a)%p;
		a = (a+a)%p;
		b >>= 1;
	}
	return ans;
}
ll qpow(ll a, ll b, ll p)
{
	ll ans = 1;
	while(b){
		if(b & 1LL)	ans = multi(ans, a, p);
		a = multi(a, a, p);
		b >>= 1;
	}
	return ans;
}

bool MR(ll n)
{
	if(n == 2)	return true;
	int s = 20, i, t = 0;
	ll u = n-1;
	while(!(u&1)){
		t++;
		u >>= 1;
	}
	while(s--){
		ll a = rand()%(n-2)+2;
		x[0] = qpow(a, u, n);
		for(i = 1; i <= t; i++){
			x[i] = multi(x[i-1], x[i-1], n);
			if(x[i] == 1 && x[i-1] != 1 && x[i-1] != n-1)	return false;	
		}
		if(x[t] != 1)	return false;
	}
	return true;	
}

ll gcd(ll a, ll b)
{
	if(b == 0)	return a;
	else	return gcd(b, a%b);
}

ll Pollard_Rho(ll n, int c)
{
	ll i = 1, k = 2, x = rand()%(n-1)+1, y = x;
	while(1){
		i++;
		x = (multi(x, x, n)+c)%n;
		ll p = gcd((y-x+n)%n, n);
		if(p != 1 && p != n)	return p;
		if(y == x)	return n;
		if(i == k){
			y = x;
			k <<= 1;
		}	
	}
}

void find(ll n, int c)
{
	if(n == 1)	return;
	if(MR(n)){
		aria.push(n);
		return;
	}
	ll p = n, k = c;
	while(p >= n){
		p = Pollard_Rho(p, c--);
	}
	find(p, k);
	find(n/p, k);
}

int main()
{
	ll n;
	while(~scanf("%lld", &n)){
		find(n, 107);
		cout << aria.front();
		aria.pop();
		while(!aria.empty()){
			cout << "*" << aria.front();
			aria.pop();
		}
		cout << endl;
	}
	return 0;
}
POJ 1811這是一道裸題,輸出最小素數即可