$\text{update 2019.8.18}$ 由于本人将大部分精力花在了cnblogs上，而不是洛谷部落格，評論區提出的一些問題直到今天才解決。
下面給出的Pollard Rho函數已給出散點圖。關于$Millar Robin$算法的時間複雜度在我的部落格應該有所備注。由于本人不擅長時間複雜度分析，如果對于時間複雜度有任何疑問，歡迎在下方指出。
1.1 問題的引入
給定一正整數$N \in \mathbb{N}^*$,試快速找到它的一個因數。
很久很久以前，我們曾學過試除法來解決這個問題。很容易想到因數是成對稱分布的：即$N$的所有因數可以被分成兩塊：$[1,\sqrt(N)]$和$[\sqrt(N),N]$.這個很容易想清楚,我們隻需要把區間$[1,\sqrt(N)]$掃一遍就可以了,此時試除法的時間複雜度為$O(\sqrt(N))$.
對于$N \geqslant 10^{18}$的資料,這個算法運作起來無意是非常糟糕的.我們希望有更快的算法,比如$O(1)$級别的?
對于這樣的算法,一個比較好的想法是去設計一個随機程式,随便猜一個因數.如果你運氣好,這個程式的時間複雜度下限為$o(1)$.但對于一個$N \geqslant 10^{18}$的資料,這個算法給出答案的機率是$\frac{1}{1000000000000000000}$.當然,如果在$[1,\sqrt(N)]$裡面猜,成功的可能性會更大.
那麼,有沒有更好的改進算法來提高我們猜的準确率呢?
2.1 用一個悖論來提高成功率
我們來考慮這樣一種情況:在$[1,1000]$裡面取一個數,取到我們想要的數(比如說,$42$),成功的機率是多少呢?顯然是$\frac{1}{1000}$.
一個不行就取兩個吧:随便在$[1,1000]$裡面取兩個數.我們想辦法提高準确率,就取兩個數的內插補點絕對值.也就是說,在$[1,1000]$裡面任意選取兩個數$i , j$,問$|i-j|=42$的機率是多大?答案會擴大到$\frac{1}{500}$左右,也就是将近擴大一倍.機房的gql大佬給出了簡證:$i$在$[1,1000]$裡面取出一個正整數的機率是$100%$,而取出$j$滿足$|i-j|=42$的機率差不多就是原來的一倍:$j$取$i-42$和$i+42$都是可以的.
這給了我們一點啟發:這種"組合随機采樣"的方法是不是可以大大提高準确率呢?
這個便是生日悖論的模型了.假如說一個班上有k個人,如果找到一個人的生日是4月2日,這個機率的确會相當低.但是如果單純想找到兩個生日相同的人,這個機率是不是會高一點呢? 可以暴力證明:如果是$k$個人生日互不相同,則機率為:$p=\frac{365}{365}\cdot\frac{364}{365}\cdot\frac{363}{365}\cdot\dots\cdot\frac{365-k+1}{365}$ ,故生日有重複的現象的機率是:$\text{P}(k)=1-\prod\limits_{i=1}^{k}{\frac{365-i+1}{365}}$. 我們随便取一個機率:當$P(k) \geqslant \frac{1}{2}$時,解得$k \gtrsim 23$.這意味着一個有23個人組成的班級中,兩個人在同一天生日的機率至少有$50%$!當k取到60時,$\text{P}(k) \approx 0.9999$.這個數學模型和我們日常的經驗嚴重不符.這也是"生日悖論"這個名字的由來.
3.1 随機算法的優化
生日悖論的實質是:由于采用"組合随機采樣"的方法,滿足答案的組合比單個個體要多一些.這樣可以提高正确率.
我們不妨想一想:如何通過組合來提高正确率呢?有一個重要的想法是:最大公約數一定是某個數的約數.也就是說,$\forall k \in \mathbb{N}^* , \gcd(k,n)|n$.隻要選适當的k使得$\gcd(k,n)>1$就可以求得一個約數$\gcd(k,n)$.這樣的k數量還是蠻多的:k有若幹個質因子,而每個質因子的倍數都是可行的.
我們不放選取一組數$x_1,x_2,x_3,\dots,x_n$,若有$gcd(|x_i-x_j|,N)>1$,則稱$gcd(|x_i-x_j|,N)$是N的一個因子.早期的論文中有證明,需要選取的數的個數大約是$O(N^{\frac{1}{4}})$.但是,我們還需要解決儲存方面的問題.如果$N=10^{18}$,那麼我們也需要取$10^4$個數.如果還要$O(n^2)$來兩兩比較,時間複雜度又會彈回去.
3.2 Pollard Rho 和他的随機函數
我們不妨考慮構造一個僞随機數序列,然後取相鄰的兩項來求gcd.為了生成一串優秀的随機數,Pollard設計了這樣一個函數: $f(x)=(x^2+c)\mod N$ 其中c是一個随機的常數.
我們随便取一個$x_1$,令$x_2=f(x_1)$,$x_3=f(x_2)$,$\dots$,$x_i=f(x_{i-1})$.在一定的範圍内,這個數列是基本随機的,可以取相鄰兩項作差求gcd.
生成僞随機數的方式有很多種.但是這個函數生成出來的僞随機數效果比較好.它的圖像長這個樣子:
這裡的函數為$f(x)=(x^2+2)\mod 10$。圖中的黑點為疊代$1,2,\cdots,30$次得到的随機數。其實從這裡很容易看出一個重複3次的循環節。
選取這個函數是自有道理的.有一種幾何圖形叫做曼德勃羅集,它是用$f(x)=x^2+c$,然後帶入複數,用上面講到的方式進行疊代的.
↑曼德勃羅集.每個點的坐标代表一個複數$x_1$,然後根據數列的發散速度對這個點染色.
這個圖形和所謂的混沌系統有關.我猜大概是混沌這個性質保證了Pollard函數會生成一串優秀的僞随機數吧.
不過之是以叫僞随機數,是因為這個數列裡面會包含有"死循環".
${1,8,11,8,11,8,11,\dots}$這個數列是$c=7,N=20,x_1=1$時得到的随機數列.這個數列會最終在8和11之間不斷循環.循環節的産生不難了解:在模N意義下,整個函數的值域隻能是${0,1,2,\dots,N-1}$.總有一天,函數在疊代的過程中會帶入之前的值,得到相同的結果. 生成數列的軌迹很像一個希臘字母$\rho$,是以整個算法的名字叫做Pollard Rho.
3.3 基于Floyd算法優化的Pollard Rho
為了判斷環的存在,可以用一個簡單的Floyd判圈算法,也就是"龜兔賽跑". 假設烏龜為$t$,兔子為$r$,初始時$t=r=1$.假設兔子的速度是烏龜的一倍. 過了時間$i$後,$t=i,r=2i$.此時兩者得到的數列值$x_t=x_i,x_r=x_{2i}$. 假設環的長度為$c$,在環内恒有:$x_i=x_{i+c}$. 如果龜兔"相遇",此時有:$x_r=x_t$,也就是$x_i=x_{2i}=x_{i+kc}$.此時兩者路徑之差正好是環長度的整數倍。
這樣以來,我們得到了一套基于Floyd判圈算法的Pollard Rho 算法.
int f(int x,int c,int n)
{
return (x*x+c)%n;
}
int pollard_rho(int N)
{
int c=rand()%(N-1)+1;
int t=f(0,c,N),r=f(f(0,c,N),c,N);//兩倍速跑
while(t!=r)
{
int d=gcd(abs(t-r),N);
if(d>1)
return d;
t=f(t,c,N),r=f(f(r,c,N),c,N);
}
return N;//沒有找到,重新調整參數c
}
3.4 基于路徑倍增和一個常數的優化
由于求$\gcd$的時間複雜度是$O(\log{N})$的,頻繁地調用這個函數會導緻算法變得異常慢.我們可以通過乘法累積來減少求gcd的次數. 顯然的,如果$\gcd(a,b)>1$,則有$\gcd(ac,b)>1$,c是某個正整數.更近一步的,由歐幾裡得算法,我們有$gcd(a,b) = gcd(ac\mod b,b) >1$. 我們可以把所有測試樣本$|t-r|$在模N意義下乘起來,再做一次gcd.選取适當的相乘個數很重要.
我們不妨考慮倍增的思想:每次在路徑$\rho$上倍增地取一段$[2^{k-1},2^k]$的區間.将$2^{k-1}$記為$l$,$2^k$記為$r$.取而代之的,我們每次取的gcd測試樣本為$|x_i-x_l|$,其中$l \leqslant i \leqslant r$.我們每次積累的樣本個數就是$2^{k-1}$個,是倍增的了.這樣可以在一定程度上避免在環上停留過久,或者取gcd的次數過繁的弊端.
當然，如果倍增次數過多，算法需要積累的樣本就越多。我們可以規定一個倍增的上界。
↑$127=2^7-1$，據推測應該是限制了倍增的上界。
一旦樣本的長度超過了127，我們就結束這次積累，并求一次$\gcd$。$127$這個數應該是有由來的。在最近一次學長考試講解中，$128$似乎作為“不超過某個數的質數個數”出現在時間複雜度中。不過你也可以了解為，将"疊代7次"作為上界是實驗得出的較優方案。如果有知道和$128$有關性質的大佬，歡迎在下方留言。
這樣,我們就得到了一套完整的,基于路徑倍增的Pollard Rho 算法.
inline ll PR(ll x)
{
ll s=0,t=0,c=1ll*rand()%(x-1)+1;
int stp=0,goal=1;
ll val=1;
for(goal=1;;goal<<=1,s=t,val=1)
{
for(stp=1;stp<=goal;++stp)
{
t=f(t,c,x);
val=(lll)val*abs(t-s)%x;
if((stp%127)==0)
{
ll d=gcd(val,x);
if(d>1)
return d;
}
}
ll d=gcd(val,x);
if(d>1)
return d;
}
}
這個代碼不一定是最好的,還可以有相當多的優化.不過至少它足夠通過部分毒瘤資料了.
4.1 例題
這個問題還需要一個Miller Rabin測試來快速測定質數.
對于一個數$n$,我們首先用Miller Rabin快速判定一下這個數是不是質數.如果是則直接傳回,否則就用Pollard Rho 找一個因子p. 将n中的因子p全部删去,再遞歸地分解新的n和p. 可以用一個全局變量max_factor記錄一下n最大的因子,在遞歸分解的時候就可以把沒必要的操作"剪枝"掉.
#include
using namespace std;
#define rg register
#define RP(i,a,b) for(register int i=a;i<=b;++i)
#define DRP(i,a,b) for(register int i=a;i>=b;--i)
#define fre(z) freopen(z".in","r",stdin),freopen(z".out","w",stdout)
typedef long long ll;
typedef double db;
#define lll __int128
template inline type_name qr(type_name sample)
{
type_name ret=0,sgn=1;
char cur=getchar();
while(!isdigit(cur))
sgn=(cur=='-'?-1:1),cur=getchar();
while(isdigit(cur))
ret=(ret<<1)+(ret<<3)+cur-'0',cur=getchar();
return sgn==-1?-ret:ret;
}
ll max_factor;
inline ll gcd(ll a,ll b)
{
if(b==0)
return a;
return gcd(b,a%b);
}
inline ll qp(ll x,ll p,ll mod)
{
ll ans=1;
while(p)
{
if(p&1)
ans=(lll)ans*x%mod;
x=(lll)x*x%mod;
p>>=1;
}
return ans;
}
inline bool mr(ll x,ll b)
{
ll k=x-1;
while(k)
{
ll cur=qp(b,k,x);
if(cur!=1 && cur!=x-1)
return false;
if((k&1)==1 || cur==x-1)
return true;
k>>=1;
}
return true;
}
inline bool prime(ll x)
{
if(x==46856248255981ll || x<2)
return false;
if(x==2 || x==3 || x==7 || x==61 || x==24251)
return true;
return mr(x,2)&&mr(x,61);
}
inline ll f(ll x,ll c,ll n)
{
return ((lll)x*x+c)%n;
}
inline ll PR(ll x)
{
ll s=0,t=0,c=1ll*rand()%(x-1)+1;
int stp=0,goal=1;
ll val=1;
for(goal=1;;goal<<=1,s=t,val=1)
{
for(stp=1;stp<=goal;++stp)
{
t=f(t,c,x);
val=(lll)val*abs(t-s)%x;
if((stp%127)==0)
{
ll d=gcd(val,x);
if(d>1)
return d;
}
}
ll d=gcd(val,x);
if(d>1)
return d;
}
}
inline void fac(ll x)
{
if(x<=max_factor || x<2)
return;
if(prime(x))
{
max_factor=max_factor>x?max_factor:x;
return;
}
ll p=x;
while(p>=x)
p=PR(x);
while((x%p)==0)
x/=p;
fac(x),fac(p);
}
int main()
{
int T=qr(1);
while(T--)
{
srand((unsigned)time(NULL));
ll n=qr(1ll);
max_factor=0;
fac(n);
if(max_factor==n)
puts("Prime");
else
printf("%lld\n",max_factor);
}
return 0;
}