素性測試

測試一個數是否為素數。

當這個數很小時，我們當然可以用試除法（用2到sqrt(N)看能否整除N）來做素性測試。

當這個數很大時，例如1024bit即1~2^1024，這個方法就顯得效率太低。考慮一個數有k位，每增加一位，N增大一倍，試除法此時是指數級別的算法。

對於大整數的素性測試，一般用Miller-Rabin算法。它是一個基於概率的算法，是費馬小定理（若n是一個素數，a^(n-1) mod n == 1；反之，n有可能是一個素數）的一個改進。

理論基礎：若n是一個素數，a<n 滿足 a^2 mod n = 1 當且僅當 a mod n = 1 或 a mod n = -1（即a mod n = n-1）。取n = q * 2^k，即若n為素數以下兩個條件必有一個滿足：1. a^q ≡ 1(mod n)；2. for j = 0~k, a^(q*2^j) ≡ n-1 (mod n)。（因為正負1的平方必為1）

更詳細的理論解釋請看CANS的P252.http://book.douban.com/subject/5372307/

另外，為了提高效率，這個算法需要用到快速冪取模算法。http://www.cnblogs.com/7hat/p/3398394.html

算法的時間復雜度比試除法顯然要高得多。注意到算法只有一個循環系數k，調用的快速冪取模也是和k有關，是多項式級別的算法。

下面給出python代碼，是因為python直接支持大整數，看起來會更加簡單。

import random

"""
e = e0*(2^0) + e1*(2^1) + e2*(2^2) + ... + en * (2^n)

b^e = b^(e0*(2^0) + e1*(2^1) + e2*(2^2) + ... + en * (2^n))
    = b^(e0*(2^0)) * b^(e1*(2^1)) * b^(e2*(2^2)) * ... * b^(en*(2^n)) 

b^e mod m = ((b^(e0*(2^0)) mod m) * (b^(e1*(2^1)) mod m) * (b^(e2*(2^2)) mod m) * ... * (b^(en*(2^n)) mod m) mod m
"""
def fastExpMod(b, e, m):
    result = 1
    while e != 0:
        if (e&1) == 1:
            # ei = 1, then mul
            result = (result * b) % m
        e >>= 1
        # b, b^2, b^4, b^8, ... , b^(2^n)
        b = (b*b) % m
    return result

def primeTest(n):
    q = n - 1
    k = 0
    #Find k, q, satisfied 2^k * q = n - 1
    while q % 2 == 0:
        k += 1;
        q /= 2
    a = random.randint(2, n-2);
    #If a^q mod n= 1, n maybe is a prime number
    if fastExpMod(a, q, n) == 1:
        return "inconclusive"
    #If there exists j satisfy a ^ ((2 ^ j) * q) mod n == n-1, n maybe is a prime number
    for j in range(0, k):
        if fastExpMod(a, (2**j)*q, n) == n - 1:
            return "inconclusive"
    #a is not a prime number
    return "composite"
print primeTest(93450983094850938450983409621);