前言

Miller_Rabin是用來快速判斷質數，而Pollard_Rho則是用來快速分解質因數
當然兩個演算法都很玄學，尤其是Pollard_Rho其實是Bogo失散多年的兄弟
Pollard_Rho很不穩定，例如當n=134579128597851時程式可以跑0.1s~1.0s不等

費馬小定理

當p為質數且a和p互質時，有
$a^{p-}$

1 ≡ 1   ( m o d   p ) a^{p-1}\equiv1\,(mod\,p) 一定成立
（當p為合數時可能成立，如果不成立就一定為合數）
現在要證明這個式子

首先由於p是質數，所以有一個優美的性質：

m o d   p ) xy\equiv 0\,(mod\,p) (x,y>0)成立時，x和y中必須要有至少一個為p的倍數
顯然
$(x\,mod\,p)(y\,mod\,p)=kp$，則
$\frac{(x\,mod\,p)(y\,mod\,p)}{p}=k$
如果x和y mod p都不為0，則找不出能湊成質數p的因子
（當然如果p是合數就不同了，例如8*9 mod 12=0）

---

考慮a,2a,3a,…,(p-1)a在mod p意義下的值
有兩個性質
①這p-1個值兩兩不同
反證法證明：
假設ma和na（m≠n）在mod p意義下相同，則
$ma\equiv na\,(mod\,p)$
則
$(m-n)a\equiv 0\,(mod\,p)$
因為ap互質，所以只可能(m-n)|p
而因為1≤m,n≤p-1且m≠n，所以m-n不可能是p的倍數

②這p-1個值中不會有0
反證法：
假設某個xa mod p=0，則
$xa\equiv 0\,(mod\,p)$
因為ap互質，所以只可能x mod p=0
因為1≤x≤p-1，所以x mod p≠0

因為各不相同且一共只有p-1個數可選，所以a~(p-1)a mod p這p-1個數只可能把1~p-1的數各取一遍
所以
$a*2a*3a*...*(p-1)a\equiv 1*2*3*...*(p-1)\,(mod\,p)$
約掉 $(p-1)!$
$a^{p-1}\equiv 1\,(mod\,p)$
得證。

Miller_Rabin

由於費馬小定理實在太好♂用了，所以自然會有人想到用其來判斷質數
於是xjb找了一堆合數p和一堆底數a來判斷後，發現某些合數也可以通過費馬小定理的測試
但是人們發現，當用來判斷的a越多時，出錯的概率就越小
於是有人猜想：如果把所有小於p的a（ap互質）都測試一遍，則判斷出來的一定是質數！

然後很快就被打臉了。
事實上，能通過所有與其互質的最小合數很小，僅為561。
（當然如果把所有小於p的a都測一遍且不考慮互質的話絕對不止這個數，而且也不會判錯，因為p與所有小於它的正整數都互質，但這樣搞還不如√n判斷）

---

所以還是考慮用某個a來判斷p
然後有一個著名的二次探測定理，也是Miller_Rabin的精髓所在。
設 $x^2\equiv1\,(mod\,p)$，則顯然
$(x-1)(x+1)\equiv0\,(mod\,p)$
若x≠1和p-1，則p不為質數

反證法證明：
若x≠1和p-1且p為質數，則(x-1)和(x+1)mod p都不為0
那麼就變成了上面的問題，ab mod p=0只有a或b中有一個為p的倍數（mod p=0）
所以(x-1)(x+1) mod p≠0
則x只能為1或p-1，與條件衝突
所以當x≠1和p-1時，p不為質數
注意這並不代表x=1或p-1時p就為質數，因為二次探測定理同樣是不確定演算法，只能判掉合數

---

這個定理有什麼用呢？
首先根據費馬小定理，a^p-1 mod p=1
然後對a^p-1不斷開方，直到有一個結果不為1或者不能繼續開方

如果結果為p-1或是當前指數為奇數，則這個數可能是質數
否則一定不是質數

---

對於每個數x，如果x=1且√x≠1，而且
(√x)² mod p=1(x)，所以
(√x-1)(√x+1) mod p=1
而√x≠1，則不滿足二次探測定理，所以p一定不為質數

但如果√x=p-1，則並不能判斷x為合數
事實上，還可以繼續開方直到某個x為1位置
但顯然，如果前面有某個x為1，則x的若干次平方必然為1
而當前的x=p-1，所以前面不可能存在某個x=1，因此繼續判斷也沒有意義了

實現

一個很顯然的方法就是把a^p-1算出來然後開方
當然這樣也很顯然meaningless

還有一個很顯然的方法就是算出a^x（x^k=p-1），然後再不斷平方直到指數為p-1，把中途的值都記錄下來，最後反著判斷
這樣做沒有問題，但是比較麻煩

事實上，可以不用記錄所有的值，只需要依次判斷當前值就行了
首先如果a^x=1，那麼顯然平方後也是1
如果當前為1，則判斷上一個數是否為p-1，不是則為合數
如果直到a^p-1都不為1，則a^(p-1)/2肯定不為p-1，不必特判

當然如果a是個質數且p=a的話要特判，因為a^a-1 mod a=0

---

Miller_Rabin的正確率很高，如果用2、3、5、7四個數作為a的話，1~25000000000000以內只有3215031751會出錯
如果輸入的數很大，則可以用慢速乘來搞

code

#include <iostream>
#include <cstdlib>
#include <cstdio>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
using namespace std;

unsigned long long n;

unsigned long long qpower(unsigned long long a,int b)
{
	unsigned long long ans=1;
	
	while (b)
	{
		if (b&1)
		ans=ans*a%n;
		
		a=a*a%n;
		b>>=1;
	}
	
	return ans;
}

bool Miller_Rabin(long long t,int N)
{
	if (t==1) return 0;
	if (t==2 || t==3 || t==5 || t==7) return 1;
	
	unsigned long long p=t-1;
	unsigned long long s,ls=0;
	
	while (!(p&1)) p>>=1;
	
	s=qpower(N,p);
	if (s==1) return 1;
	
	while (p<=t-1 && s!=1)
	{
		ls=s;
		s=s*s%t;
		p<<=1;
	}
	
	return ls==t-1;
}

bool pd(long long n)
{
	if ((n/100000==32150) && (n%100000==31751))
	return 0;
	else
	if (Miller_Rabin(n,2) && Miller_Rabin(n,3) && Miller_Rabin(n,5) && Miller_Rabin(n,7))
	return 1;
	else
	return 0;
}

int main()
{
	scanf("%lld",&n);
	
	if (pd(n))
	printf("Yes\n");
	else
	printf("No\n");
	
	//in the scope of 1~25000000000000,only 3215031751 is wrong
}

Pollard_Rho

考慮最原始分解質因數的方法：
每次列舉一個≤√n的數，將其全部除去，若最後剩下的數>1的則為其一個質因子
然而這樣的複雜度是穩定O(√n)的，考慮更快的方法

---

先考慮如何求n的一個因子
有一個很顯然的方法：
每次隨機random一個2~n-1的數，判斷是否為n的因子，是則輸出，不是則繼續
若n=p*q（pq均為質數），則期望時間複雜度為O(n)
而且c++用的是偽隨機數，有可能陷入迴圈無法自拔

生日悖論

有一個問題：
給出n個人，求n個人當中有任意兩個人生日相同的期望
證明不會，結論是當n>√365時，相同的概率>50%
可以推廣一下，從[1~x]中選n個數，當n>√x時概率>50%

對於生日悖論來說，選擇的兩個數的關係是a=b
而對於b-a=p這種情況，實際就相當於a加了某個數，結論還是一樣
所以可以推廣：從[1~x]中選n個數，當n>√x時b-a=p的概率>50%
而n=p*q（最壞情況），則只需要找√n個數然後兩兩做差，判斷差是否為n的因數，這樣成功的概率>50%

然而——
這樣做的時間複雜度實際上是O(√n²)=O(n)，而且正確率堪憂

真·Pollard_Rho

因為找b-a=p的概率實在是太小了，只有 $\frac{1}{n}$
但如果我們改一下：
判斷gcd(b-a,n)>1呢？（b>a）
顯然如果大於1，那麼gcd肯定是n的一個因數（最大公因數也是其中一個的因數）
那麼可以發現，原來b-a只可能為p或q
而現在b-a可以為p,2p,3p,…(q-1)p,q,2q,3q,…(p-1)q，因為這些數都包含了p或q
所以一共有p+q-2個數，也就是概率為 $\frac{p+q-2}{n}$
因為n=p*q，所以顯然pq越接近它們的和就越小
最壞情況p≈√n，則概率為 $\frac{√n+√n-2}{n}$，也就是 $\frac{√n}{n}$

一下提升了√n倍，所以顯然只需要找 $\sqrt[4]{n}$個數（原來是 $\sqrt{n}$個），所以兩兩匹配的複雜度為 $O(\sqrt{n}*log\,n)$，log n太小了所以一般不考慮

具體實現

但是還有一個問題：
如果n特別大，連 $\sqrt[4]{n}$都存不下要怎麼搞？
顯然可以發現，如果把數全部存下來再搞，就相當於每次隨機兩個數然後做差判斷
於是可以每次只存兩個數，不斷做差判斷即可
這樣和上面是類似的，最壞時間複雜度為 $O(\sqrt{n})$

對於隨機數的選擇，有一個很好的函式f(x)=(x²+a)%p（a可以為任意數，也可以隨機）
顯然這樣搞實在模p意義下，難免會出現環的情況
於是就有一個神奇的Floyd判環法：
設ab兩個初值相同的數，讓b按照a的兩倍速度來跑，如果在之後的某個時刻a=b則出現了環
然而這樣在某些奇♂妙的數（例如4的倍數）下會陷入死迴圈，所以可以先搞掉質因數2357再說
也可以每次都隨機一個a

然而這樣搞還是跟Miller_Rabin沒有什麼關係

---

假設給出的數是一個大質數，則Pollard_Rho可能會掛掉
所以就需要Miller_Rabin判斷當前的數是否是質數
如果要把一個數分解質因數，那麼就遞迴分解到不是質數位置

code

#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <ctime>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define abs(x) ((x>0)?(x):-(x))
using namespace std;

long long ans[101];
unsigned long long P,p,A;
int len,i;

unsigned long long chen(unsigned long long a,int b,long long p)
{
	unsigned long long ans=0;
	
	while (b)
	{
		if (b&