hihocoder 數論一·Miller-Rabin質數測試

題目1 : 數論一·Miller-Rabin質數測試

時間限制: 10000ms

單點時限: 1000ms

記憶體限制: 256MB

描述

小Hi和小Ho最近突然對密碼學産生了興趣，其中有個叫RSA的公鑰密碼算法。RSA算法的計算過程中，需要找一些很大的質數。

小Ho：要如何來找出足夠大的質數呢？

小Hi：我倒是有一個想法，我們可以先随機一個特别大的初始奇數，然後檢查它是不是質數，如果不是就找比它大2的數，一直重複，直到找到一個質數為止。

小Ho：這樣好像可行，那我就這麼辦吧。

過了一會兒，小Ho拿來了一張寫滿數字的紙條。

小Ho：我用程式随機生成了一些初始數字，但是要求解它們是不是質數太花時間了。

小Hi：你是怎麼做的啊？

說着小Hi接過了小Ho的紙條。

小Ho：比如說我要檢測數字n是不是質數吧，我就從2開始枚舉，一直到sqrt(n)，看能否被n整除。

小Hi：那就對了。你看紙條上很多數字都是在15、16位左右，就算開方之後，也有7、8位的數字。對于這樣大一個數字的循環，顯然會很花費時間。

小Ho：那有什麼更快速的方法麼？

小Hi：當然有了，有一種叫做Miller-Rabin質數測試的算法，可以很快的判定一個大數是否是質數。

​​提示：Miller-Rabin質數測試​​  

輸入

第1行：1個正整數t，表示數字的個數，10≤t≤50

第2..t+1行：每行1個正整數，第i+1行表示正整數a[i]，2≤a[i]≤10^18

輸出

第1..t行：每行1個字元串，若a[i]為質數，第i行輸出"Yes"，否則輸出"No"

樣例輸入

3
3
7
9      

樣例輸出

Yes
Yes
No      

提示：Miller-Rabin質數測試

小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：好！這樣就能夠順利的找到大質數了。

本題的提示參考了Matrix67的Blog​和wikipedia的詞條。

Matrix67的Blog有更多的細節描寫。Wiki中的僞代碼比上文中的簡潔一些，并且有介紹了一些小技巧：比如如果n<2^64，隻用選取a=2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37做測試即可

AC代碼：

#include<bits/stdc++.h>
#pragma comment(linker, "/STACK:1024000000,1024000000")


using namespace std;
#define eps 1e-9
#define PI acos(-1.0)
#define INF 0x3f3f3f3f
#define LLINF 1LL<<50
#define speed std::ios::sync_with_stdio(false);

typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;
typedef pair<ll, ll> pll;
typedef complex<ld> point;
typedef pair<int, int> pii;
typedef pair<pii, int> piii;
typedef vector<int> vi;

#define CLR(x,y) memset(x,y,sizeof(x))
#define CPY(x,y) memcpy(x,y,sizeof(x))
#define clr(a,x,size) memset(a,x,sizeof(a[0])*(size))
#define cpy(a,x,size) memcpy(a,x,sizeof(a[0])*(size))
#define debug(a) cout << #a" = " << (a) << endl;
#define debugarry(a, n) for (int i = 0; i < (n); i++) { cout << #a"[" << i << "] = " << (a)[i] << endl; }

#define mp(x,y) make_pair(x,y)
#define pb(x) push_back(x)
#define lowbit(x) (x&(-x))

#define MID(x,y) (x+((y-x)>>1))
#define getidx(l,r) (l+r | l!=r)
#define ls getidx(l,mid)
#define rs getidx(mid+1,r)
#define lson l,mid
#define rson mid+1,r

template<class T>
inline bool read(T &n)
{
    T x = 0, tmp = 1;
    char c = getchar();
    while((c < '0' || c > '9') && c != '-' && c != EOF) c = getchar();
    if(c == EOF) return false;
    if(c == '-') c = getchar(), tmp = -1;
    while(c >= '0' && c <= '9') x *= 10, x += (c - '0'),c = getchar();
    n = x*tmp;
    return true;
}
template <class T>
inline void write(T n)
{
    if(n < 0)
    {
        putchar('-');
        n = -n;
    }
    int len = 0,data[20];
    while(n)
    {
        data[len++] = n%10;
        n /= 10;
    }
    if(!len) data[len++] = 0;
    while(len--) putchar(data[len]+48);
}

ll prime[6] = {2, 3, 5, 233, 331};
ll qmul(ll x, ll y, ll mod)   // 乘法防止溢出， 如果p * p不爆LL的話可以直接乘； O(1)乘法或者轉化成二進制加法(快速加)
{
    ll ret = 0;
    while(y) {
        if(y & 1)
            ret = (ret + x) % mod;
        x = x * 2 % mod;
        y >>= 1;
    }
    return ret;
}
ll qpow(ll a, ll n, ll mod)
{
    ll ret = 1;
    while(n)
    {
        if(n & 1) ret = qmul(ret, a, mod);
        a = qmul(a, a, mod);
        n >>= 1;
    }
    return ret;
}
bool Miller_Rabin(ll p)
{
    if(p < 2) return 0;
    if(p != 2 && p % 2 == 0) return 0;
    ll s = p - 1;
    while(! (s & 1)) s >>= 1;
    for(int i = 0; i < 5; ++i)
    {
        if(p == prime[i]) return 1;
        ll t = s, m = qpow(prime[i], s, p);
        while(t != p - 1 && m != 1 && m != p - 1)
        {
            m = qmul(m, m, p);
            t <<= 1;
        }
        if(m != p - 1 && !(t & 1)) return 0;
    }
    return 1;
}

int main()
{
    int n;
    scanf("%d",&n);
    while(n--)
    {
        ll num;
        scanf("%lld",&num);
        if(Miller_Rabin(num))
            puts("Yes");
        else
            puts("No");
    }
    return 0;
}