演算法記錄——pallard rho（大整數分解）
要解決的問題很簡單，對一個整數進行分解質因數。
首先還是效率非常低的暴力演算法，相信大家都會，不多提。
和上次一樣，當數達到非常大的時候，分解將變得非常困難。於是這次又帶來一個提升分解速度的“非完美”演算法。之所以打引號，是因為這次不完美的不是結果，而是時間效率。
Pollard Rho演算法分解一個數n的過程大體上是這樣子的：
1、找到一個數p，使得p|n
，將n分解為p與n/p
2、如果p或n/p不為質數，將其帶入遞迴上述過程
3、如果其是質數，將其記錄並退出
是不是很sb。。。有人就會問了：這跟暴力分解有什麼區別？好像時間複雜度還比暴力高一些。。。
所以：下面的優化才是關鍵。
第一個優化，使用Miller Rabin演算法判定其是否為質數，這個不多提。
關鍵就在於接下來的這個優化。
對於一個大整數n，我們要找到一個p滿足p|n
，這顯然如大海撈針。但是如果我們要找出p1、p2，使得abs(p1−p2)|n
,這看起來似乎要容易一些。實際上我們只需要找出gcd(abs(p1−p2),n)>1
的p1、p2，則其gcd值肯定為n的約數。這看起來又容易了一些。
實際上，不止容易一些，而是容易許多。根據某個玄學理論（生日悖論，詳見百度，在此不贅述），這種兩兩比較的方式，在加入比較的數越來越多的時候，其概率會大大提升，比找一個數的概率提升快很多。
於是現在，找p的過程變成了這個樣子：
1、找到一個數p1
2、通過某種玄學推導手段找出一個與p1對應的p2
3、判斷gcd(abs(p1−p2),n)
是否大於1，不大於則將p2作為新的p1，重複過程，否則就找到了
怎麼又是玄學?因為只有通過推導手段，才能保證不做重複判斷，浪費時間。理論上的推導手段可以有很多，但實際使用中統一使用如下公式推導：
p2=(p12+c)modn
其中c為隨機常數。
這個公式的好處:
1、顯然推匯出來的p2-p1差值基本不會相等。
2、可以證明，該推導結果會出現迴圈。也就是說，在出現迴圈之前，結果不會重複，少做了許多無用的判斷。
出現迴圈了怎麼辦？換一個隨機常數再搞。這就是該演算法“非完美”的地方，如果人品太差那就。。。不過根據上面函式影象可知，兩個隨機常數產生的推導結果基本不會有重複，所以就可以放心開搞了。
最後一點，判環怎麼判？floyd判圈演算法搞定。（一個標記以另一個標記幾倍速度走，在環上總能碰到。詳見百度）
需要注意的是，之所以不能一個標記定在原地，是因為迴圈節不一定在開頭就產生，可能走著走著才遇到迴圈。這條路徑就類似於ρ
，Pollard Rho演算法因此得名。
最後附上程式碼。
C++示例程式碼
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <stdio.h>
#include <algorithm>
#include <math.h>
#include <stdlib.h>
#include<time.h>
#define ll long long
#define INF 0x3f3f3f3f
#define maxn 10000+10
#define cle(a) memset(a,0,sizeof(a))
const double eps=1e-5;
using namespace std;
const int S=20;//隨機演算法判定次數，S越大，判錯概率越小
//計算 (a*b)%c. a,b都是long long的數，直接相乘可能溢位的
// a,b,c <2^63
ll mult_mod(ll a,ll b,ll c)
{
a%=c;
b%=c;
ll ret=0;
while(b)
{
if(b&1)
{
ret+=a;
ret%=c;
}
a<<=1;
if(a>=c)
a%=c;
b>>=1;
}
return ret;
}
ll pow_mod(ll x,ll n,ll mod)
{
if(n==1)
return x%mod;
x%=mod;
ll tmp=x;
ll ret=1;
while(n)
{
if(n&1)
ret=mult_mod(ret,tmp,mod);
tmp=mult_mod(tmp,tmp,mod);
n>>=1;
}
return ret;
}
//以a為基,n-1=x*2^t a^(n-1)=1(mod n) 驗證n是不是合數
//一定是合數返回true,不一定返回false
bool check(ll a,ll n,ll x,ll t)
{
ll ret=pow_mod(a,x,n);
ll last=ret;
for(int i=1; i<=t; i++)
{
ret=mult_mod(ret,ret,n);
if(ret==1&&last!=1&&last!=n-1)
return true;//合數
last=ret;
}
if(ret!=1)
return true;
return false;
}
// Miller_Rabin()演算法素數判定
//是素數返回true.(可能是偽素數，但概率極小)
//合數返回false;
bool Miller_Rabin(ll n)
{
if(n<2)
return false;
if(n==2)
return true;
if((n&1)==0)
return false;//偶數
ll x=n-1;
ll t=0;
while((x&1)==0)
{
x>>=1;
t++;
}
for(int i=0; i<S; i++)
{
ll a=rand()%(n-1)+1;//rand()需要stdlib.h標頭檔案
if(check(a,n,x,t))
return false;//合數
}
return true;
}
ll factor[100];//質因數分解結果（剛返回時是無序的）
int tol;//質因數的個數。陣列小標從0開始
ll gcd(ll a,ll b)
{
if(a==0)
return 1;//???????
if(a<0)
return gcd(-a,b);
while(b)
{
ll t=a%b;
a=b;
b=t;
}
return a;
}
ll Pollard_rho(ll x,ll c)
{
ll i=1,k=2;
ll x0=rand()%x;
ll y=x0;
while(1)
{
i++;
x0=(mult_mod(x0,x0,x)+c)%x;
ll d=gcd(y-x0,x);
if(d!=1&&d!=x)
return d;
if(y==x0)
return x;
if(i==k)
{
y=x0;
k+=k;
}
}
}
//對n進行素因子分解
void findfac(ll n)
{
if(Miller_Rabin(n))//素數
{
factor[tol++]=n;
return;
}
ll p=n;
while(p>=n)
p=Pollard_rho(p,rand()%(n-1)+1);
findfac(p);
findfac(n/p);
}
int main()
{
ll n;
while(cin>>n)
{
tol=0;
findfac(n);
for(int i=0; i<tol; i++)
cout<<factor[i]<<endl;
//質因子
}
return 0;
}
pytho示例程式碼
from random import randint
f = []
def gcd(m, n):
if n > 0:
return gcd(n, m % n)
return m
def isPrime(n):
if n < 2:
return True
for i in range(2, int(n ** 0.5) + 1):
if n % i == 0:
return False
return True
def ifenough(num):
m = 1
for i in f:
m *= i
if m == num:
return 0
if m > num:
return 1
return 2
def PollardRho(num, n):
if (isPrime(n)):
f.append(n)
return;
x = []
i = 1
x.append(-1)
x.append(randint(1, n))
y = x[1]
k = 2
while True:
i += 1
x.append((abs(x[i - 1] * x[i - 1] - 1)) % n)
d = gcd(abs(y - x[i]), n)
if d > 1 and d < n:
PollardRho(num, d)
PollardRho(num, n / d)
if i == k:
y = x[i]
k = 2 * k
if x.index(x[i], 0, i + 1) != i or ifenough(num) == 0 or ifenough(num) == 1:
break
def repeat(n):
while True:
f[:] = []
PollardRho(n, n)
if ifenough(n) == 0:
f.sort(key=None, reverse=False)
for index, j in enumerate(f):
if index != len(f) - 1:
print(j, '*')
else:
print(j)
break
if __name__ == '__main__':
n = int( input('please enter the number you want to divide:'))
repeat(n)