素性测试

测试一个数是否为素数。

当这个数很小时，我们当然可以用试除法（用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);