http://hihocoder.com/contest/hiho92/problem/1

小Hi：這種質數演算法是基於費馬小定理的一個擴充套件，首先我們要知道什麼是費馬小定理：
費馬小定理：對於質數p和任意整數a，有a^p ≡ a(mod p)(同餘)。反之，若滿足a^p ≡ a(mod p)，p也有很大概率為質數。
將兩邊同時約去一個a，則有a^(p-1) ≡ 1(mod p)
也即是說：假設我們要測試n是否為質數。我們可以隨機選取一個數a，然後計算a^(n-1) mod n，如果結果不為1，我們可以100%斷定n不是質數。
否則我們再隨機選取一個新的數a進行測試。如此反覆多次，如果每次結果都是1，我們就假定n是質數。
該測試被稱為Fermat測試。需要注意的是：Fermat測試不一定是準確的，有可能出現把合數誤判為質數的情況。
Miller和Rabin在Fermat測試上，建立了Miller-Rabin質數測試演算法。

與Fermat測試相比，增加了一個二次探測定理：
如果p是奇素數，則 x^2 ≡ 1(mod p)的解為 x ≡ 1 或 x ≡ p - 1(mod p)
如果a^(n-1) ≡ 1 (mod n)成立，Miller-Rabin演算法不是立即找另一個a進行測試，而是看n-1是不是偶數。如果n-1是偶數，另u=(n-1)/2，並檢查是否滿足二次探測定理即a^u ≡ 1 或 a^u ≡ n - 1(mod n)。
舉個Matrix67 Blog上的例子，假設n=341，我們選取的a=2。則第一次測試時，2^340 mod 341=1。由於340是偶數，因此我們檢查2^{170，得到2}170 mod 341=1，滿足二次探測定理。同時由於170還是偶數，因此我們進一步檢查2^85 mod 341=32。此時不滿足二次探測定理，因此可以判定341不為質數。
將這兩條定理合起來，也就是最常見的Miller-Rabin測試。
但一次MR測試仍然有一定的錯誤率。為了使我們的結果儘可能的正確，我們需要進行多次MR測試，這樣可以把錯誤率降低。
寫成虛擬碼為：

Miller-Rabin(n):
	If (n <= 2) Then
		If (n == 2) Then
			Return True
		End If
		Return False
	End If
	
	If (n mod 2 == 0) Then
		// n為非2的偶數，直接返回合數
		Return False
	End If
	
	// 我們先找到的最小的a^u，再逐步擴大到a^(n-1)
	
	u = n - 1; // u表示指數
	while (u % 2 == 0) 
		u = u / 2
	End While // 提取因子2
	
	For i = 1 .. S // S為設定的測試次數
 

		a = rand_Number(2, n - 1) // 隨機獲取一個2~n-1的數a
		x = a^u % n
		While (u < n) 
			// 依次次檢查每一個相鄰的 a^u, a^2u, a^4u, ... a^(2^k*u)是否滿足二次探測定理
			y = x^2 % n 
			If (y == 1 and x != 1 and x != n - 1)	// 二次探測定理
				// 若y = x^2 ≡ 1(mod n)
				// 但是 x != 1 且 x != n-1
				Return False
			End If
			x = y
			u = u * 2 
		End While
		If (x != 1) Then	// Fermat測試
			Return False
		End If
	End For
	Return True

值得一提的是，Miller-Rabin每次測試失誤的概率是1/4；進行S次後，失誤的概率是4^(-S)。
小Hi：那麼小Ho，你能計算出這個演算法的時間複雜度麼？
小Ho：恩，每一次單獨的MR測試，需要O(log n)的時間。一共要進行S次MR測試，也就是O(Slog n)。
小Hi：沒錯，這樣就能夠在很短的時間內完成質數的測試了。當然如果你還是不放心，可以把S的值設定的更高一點。
小Ho：好！這樣就能夠順利的找到大質數了。

AC程式碼：

#include <bits/stdc++.h>

using namespace std;
#define LL long long
const int Mod = 1e9 + 7;
const int maxn = 1e5 + 5;
const double eps = 0.00000001;
const int INF = 0x3f3f3f3f;
#define random(a, b) (((double)rand()/RAND_MAX)*(b-a)+a)

LL Quick_Mul(LL a, LL b, LL n) {
    LL res = 0;
    while(b) {
        if(b & 1) res = (res + a) % n;
         a = (a + a) % n;
        b >>= 1;

    }
    return res;
}

LL QuickPow(LL a, LL b, LL n) {
    LL ans = 1;
    while(b) {
        if(b & 1) ans = Quick_Mul(ans, a, n);
        a = Quick_Mul(a, a, n);
        b >>= 1;
    }
    return ans;
}

bool Miller_Rabin(LL n) {
    if(n <= 2) {
        if(n == 2) return true;
        return false;
    }
    if(n % 2 == 0) return false;
    LL u = n - 1;
    while(u % 2 == 0) u /= 2;
    int S = 100;
    srand((LL)time(0));
    for (int i = 1; i <= S; i ++){
        LL a = rand() % (n - 2)+ 2;
        LL x = QuickPow(a, u, n);
        while(u < n) {
            LL y = QuickPow(x, 2, n);
            if(y == 1 && x != 1 && x != n - 1)
                return false;
            x = y;
            u = u * 2;
        }
        if(x != 1) return false;
    }
    return true;
}

int main()
{
    int T;
    cin >> T;
    LL a[60];
    for (int i = 0; i < T; i ++) cin >> a[i];
    for (int i = 0; i < T; i ++) {
        if(Miller_Rabin(a[i])) cout << "Yes\n";
        else cout << "No\n";
    }
    return 0;
}