2017-07-19 08:54 Amphetamine:能發一下代碼嗎？
應我那位謎一樣好友的邀約，我打算好好看一看Miller-Rabin和Pollard-Rho算法。很奇怪，各種地方有很多代碼描述詳細過程，但我仍舊很懵。也許是我太弱了，不能從那些“魚龍混雜”的代碼中找出本質上的共性。那么，我們現在來討論一下吧。
首先，大整數分解現在仍然是個世界級的難題，在“公共密鑰”的研究上有着重要的作用。
！！先判斷質數！！
試除法：原始的根號算法
額。不想說了。正經一點。
Miller-Rabin：判斷一個“p”是否是質數
首先，快速冪是極為必要的。
YYR說，這種東西是要靠RP的。但是，99.999%應該還是足以解決一些問題（或是大多數）了。
a^(p-1)≡1(mod p)，gcd(a,p)=1
費馬小定理是極強的質數性質。只要找到一個a，與那個“p”不滿足上述條件，就可以跳合數了？答案是肯定的。
但是，如果找不到呢？就是質數？呵呵，很可惜啊。因為存在Carmichael數：無論取多少個a，有一個不滿足，算我輸。
比如：561 = 11*51就是一個Carmichael數。
為了99.999%，YYR給了一種相對靠譜的方法。
【PART 1：必要的第一次粗篩】
先把2,3,5,7,11,13拿來篩一遍。
這樣理論上可以篩掉：1-(1-1/2)*(1-1/3)*(1-1/5)*(1-1/7)*(1-1/11)*(1-1/13)=1-19.18%=80.82%的數。
當然，如果再加上17,19,23,29，理論上可以只會留下15.79%的數，也就是篩去了84.21%的數。
肯定，我們會發現，這個篩率會趨於收斂。
但是，讓我們專門卡素數，怎么可能只用這種方法！！！
【PART 2：費馬小定理與“二次探測定理”】
在進入這一步之前，我們還不能確定這個“p”的身份。
我們此時會先把p-1搞一搞，不停地除以2。直到發現p-1=2^k*t，而t是奇數。
這樣有何用意？且讓我們分析。有一個定理：
若a^2≡1(mod p)，則a≡±1(mod p)
YYR說到這里時，我說這很顯然，他瞪大了眼睛。因為若p|(a^2-1)，就有p|(a+1)(a-1)。這是經典的歐幾里得引理，偉大的歐幾里得就是以此證明了唯一分解定理。我們有(a+1)%p==0或(a-1)%p==0。即a≡±1(mod p)。
然后，我們可以去隨機地選取一些a。首先用快速冪算出a^t%p，然后把這個數自乘。得到了a^2t%p，然后再自乘。得到了a^4t%p，然后再自乘……我們可以發現，a^t%p自乘k次以后，就可以得到a^(2^k*t)%p，最后就是a^(p-1)%p了。
在這樣的一條路上，我們可以充分驗證上述“定理”。前面k-1個數中只要出現了1，而它的“父親”（它的“父親”自乘得到了它）不是±1，那么就可以跳合數了。最后，第k個數即a^(p-1)%p若不為1，也可以跳合數。
為什么這樣效率很高？事實證明，若“p”通過一次測試，不是素數的概率為25%。若“p”通過x次測試，則不是素數的概率為1/(4^x)。當x=5時，素數誤判的概率為1/(4^5)=1/1024=0.10%，已經直逼99.90%了。
而假如選取了4個a為2,3,5,7，則在2.5*10^13以內唯一一個失誤的數為3215031751。
所以說，選幾個a比較好？實踐證明，這種東西可以悠着來。
【小結】
1.讀入一個"p"。
2.用2,3,5,7,11,13粗篩，發現整除直接跳合數。
3.把p-1搞一搞，不停地除以2。得到奇數t，和一個數k。p-1=2^k*t。
4.隨機地找一個數a，算出x=a^t。
5.last=x,x=x*x%p,若x==1且last!=1且last!=p-1，則跳合數。該步驟重復k次。
6.現在x=p-1，若x!=1，則跳合數。
7.按心情回到第4步，99.99%與99.9999%在神犇眼里有着本質區別。
！！然后再
亂搞
！！
FERMAT：平方差
如果是偶數，就把它釜底抽薪，變成奇數。
對於一個奇數N，若不為質數，則設N=c*d，明顯c、d均為奇數。
不妨設c>d，令a=(c+d)/2,b=(c-d)/2，很明顯，N=a*a-b*b。這就是了。
我們枚舉一個a，有a*a-N=b*b，a>=sqrt(N),若a*a-N是一個完全平方數，就可以遞歸了。
但是，很明顯。這樣也是枚舉，效率很低，得不償失。
單純的POLLARD RHO：生日悖論和判圈
精彩現在繼續。有一篇文章很好：點擊這里
單純的Pollard Rho算法非常有趣，容易理解。雖然不是目前最快的算法，但它要比試除法快上多個量級。實現它的思想同樣可以用於其他地方。
該算法最早於1975年由John M. Pollard提出，而Richard Brent於1980年提出了改進版本。Wikipedia上說得特別清楚。
上一個FERMAT算法中，我們試圖分解N=c*d。那么，現在我們會繼續這一思路。假設N只有c和d兩個互異質因子，然后在[2,N]中隨機出一個數，概率是2/（N-1），很小。
但是，隨機選取23個人，他們中存在相同生日的概率大於50%。
這就是生日悖論。
但為什么是對的？如果是對的，然后呢？
【生日悖論的正確性】
我們設現在隨機選取n個人(n<=365)，希望能夠算出“存在相同生日”的概率。請注意，這個事件是指，我們能夠從中找出至少兩個生日一樣的人。
設該概率為P，很顯然它與“不存在相同生日”的概率Q之和為1。而對於Q，我們可以發現，如果人一個一個加入。對於第2個人，要讓他與第1個人不同，概率為1-1/365；對於第3個人，要讓他與第1個人和第2個人不同，概率為1-2/365；對於第4個人，要讓他與第1個人、第2個人和第3個人不同，概率為1-3/365；對於第5個人，要讓他與第1個到第4個人不同，概率為1-4/365；……；最后，對於第n個人，要讓他與第1個到第n-1個人不同，概率為1-(n-1)/365。最后，即可得Q=(1-1/365)*(1-2/365)*(1-3/365)……(1-(n-1)/365)。
中間有一個不等式：1+x<=ex(x≠0)，我不想證了。
然后，就有Q=e(-1/365)*e(-2/365)*……*e(-(n-1)/365)=e(-n(n-1)/(2*365))
於是
我們就可以算出，當n>=23時，他們中存在相同生日的概率大於50%。
其實，如果說一年有N天，那么只要n>=sqrt(N)，他們中存在相同生日的概率就會大於50%。（詳見《算導》黑書）
【然后呢？】
現在，我們隨機從[2,N]中選出n個數，當好刷出c或d的概率很小，但是如果它們兩兩作差呢？
使用類似的思路，我們可以猜（當然也可以用完全相同的證法），知道n~sqrt(N)時，概率約為50%。其實仔細想想也很有道理，因為sqrt(N)個數任意組隊，有N/2種方案。
雖然很有道理，但是不可能直接兩兩作差，因為概率還是很低。但是，你可能已經想到了。
如果是gcd(|a-b|,N)呢？於是，事情變得很好看了。
因為N=c*d，c和d都很稀有，但是c和d的倍數便是一波大軍。
所以，一個簡單的策略如下：
- 在區間[2,N-1]中隨即選取n個數，x1,x2, … … , xn
- 判斷是否存在gcd(|xi-xj| ,N) >1, 若存在，gcd(|xi-xj| ,N) 是N的一個因子 (c 或 d)
這是很可靠的。如果N有多個因子，也能用類似於FERMAT的方法遞歸解決。
但是我們算一算。
但是很早就出現了一個問題，我們需要選取大約N1/4個數，這個數量太大了，以至於我們並不能將其存放在內存中。
然后，怎么去造這n個數呢？
PollardRho說，很簡單。
【流水線】
我們發現，n個數如果先造出來，再去判斷實在太傻。於是PollardRho告訴我們，2個數就夠了。
我們可以用不斷地生成偽隨機數。然后像流水線一樣判斷。
有一個隨機函數特別棒，在此舉例：f(x)=(x^2+a)%N
這個a已經被你給定了，可以為1當然也能是rand()，只是中途一般不作變化。
代碼如下：
1 int find_factorplus(int N) { 2 a = 2; 3 for（ int i= 1; i <= 1000000; i++ ) { 4 b = f(a); 5 p = GCD( abs( b - a ) , N); 6 if( p > 1 ) return p;//Found factor: p 7 a = b; 8 } 9 return 0;//Failed. :-( 10 }
似乎很玄學，但是實際效果確實很棒。但不好的是，偽隨機數有着玄學般的循環節。
【floyd來判圈】
對於循環節。你一定會覺得，用個vis數組就好了。似乎很可以，只要你存得下，不揪心就好。
劉汝佳先生說：想象一下，假設有兩個小孩子在一個“可以無限向前跑”的跑道上賽跑，同時出發。但其中一個小孩的速度是另一個的兩倍。如果跑道是直的，跑得快的小孩永遠在前面；但如果跑道有環，則跑得快的小孩將“追上”跑得慢的小孩。
於是就這樣了。這叫做“floyd判圈”算法。
1 int find_factorplus(int N) { 2 a = 2; 3 b = a; 4 do { 5 a = f(a);//a runs once 6 b = f(f(b));//b runs twice as fast 7 p = GCD( abs( b - a ) , N); 8 if( p > 1 ) return p;//Found factor: p 9 } while( b != a ); 10 return 0;//Failed. :-( 11 }
這樣，我們就可以把退出條件溫和化，只要發現有環，那就只有退出了。而不是暴力地把i從1 for 到 1,000,000。
如果算法失敗了，我們只需要找到一個新的偽隨機函數f(x)或是一個新的a就好了。不過請放心，大多數時候你並不會失敗。
最后說一下，代碼中a的初值是2，在實際生活中，你並不需要那么講究，rand()一個也是不錯的選擇。
“最后”的POLLARD RHO：當與Miller-Rabin發生反應
我們可以發現pollard rho直到現在都還沒有與Miller-Rabin有任何聯系，但馬上就不是了。
對於pollard rho，它可以在Θ(sqrt(p))的時間復雜度內找到N的一個小因子p，這一點我們曾論證過。可見，如果N的因子很多、因子值很小的整數N來說，效率是很優異的。
但是，如果反過來呢？如果說N是大整數，恰好因子很少、因子值很大？
例如，N=2*p，p為質數。你立馬發現，N有一個因子2，然后你試圖去解決p。然后，這個很優秀的算法成了根號算法。而且直到最后，你都很難判斷這個p是否真的不可約。
但是，一旦擁有Miller-Rabin，一切便都已解決。
我們現在可以分析一下復雜度。N的質因子中，超過sqrt(N)的有且僅有一個。這樣，即使運氣極差，也能有相當的保障。
！！最后總結一下！！
斯堪福說，總結是好習慣。
對於Miller Rabin，我們需要一個快速冪，一個快速乘。先用2,3,5,7,11,13粗篩一遍，再將p的2抽盡，然后隨機地選取一些數進行二次探測與費馬小定理檢驗。
對於Pollard Rho，我們需要一個偽隨機函數f，一個常數a，一個gcd，一個abs。使用floyd判圈算法。找到一個因子后遞歸解決，中間判斷是否是質數。如果是，做記錄。
當我們在做大數質因子分解時，質因子記錄完畢后，我們常常會發現這是無序的。這就需要進行一下排序，然后離散化處理出每個質因子出現的次數。這樣就解決了。就真的解決了。
致謝
yy儒學長（NOI AG THU）
Amphetamine
IDY002
BIGBALLON
參考
劉汝佳《算法競賽入門經典》《算法競賽入門經典訓練指南》
林厚從、王宏《信息學奧賽值數學一本通》
丁堯堯idy002《信息學競賽中的數學知識小結
》
BIGBALLON《A Quick Tutorial on Pollard's Rho Algorithm》