素性测试

判断一个数$n$是否为素数有很多做法，最常见的是枚举$i$从$2$到$\lfloor \sqrt{n} \rfloor$，判断$n$是否都不能被$i$整除，代码如下：

bool isPrime(long long p){
    if(p==1)return 0;
    for(long long i=2;i*i<=p;++i)
        if(p%i==0)return 0;
    return 1;
}

然而上述算法的复杂度为$O(\sqrt{n})$，对于大数来说，这个时间是无法接受的.

为了解决大数的素性判断，有了下面的RP算法：Fermat测试.

根据Fermat小定理，若$p$为素数，则对任意$a$必有$a^{p-1} \equiv 1(mod p)$.

故对于一个$a$，若$a^{p-1} \not\equiv 1(mod p)$，则$p$为合数；若$a^{p-1} \equiv 1(mod p)$，则$p$有可能为素数，实际上当$p$为奇数时，$p$是合数的概率小于$\frac{1}{2}$.

于是对于一个奇数$p(p\geqslant 3)$和安全参数$k$：

• 随机取一个数$n(2 \leqslant n \leqslant p-1)$
• 若$(n,p)\neq 1$，则$p$为合数
• 若$n^{p-1} \not\equiv 1(mod p)$则$p$为合数，否则假定$p$为素数
• 重复上述过程$k$次

假定大数相乘的复杂度为$O(lgn \times lglgn)$，故算法复杂度为$O(k \times lg^2n \times lglgn)$，若测试得到$p$为素数，则准确率为$1-\frac{1}{2^k}$.代码如下（为了方便高精度运算，用python编写）：

import random


def GCD(a, b):
    while b:
        a, b = b, a % b
    return a


def PowMod(a, n, p):
    r, t = 1, a
    while n:
        if n & 1:
            r = (r * t) % p
        t = (t * t) % p
        n >>= 1
    return r


def FermatTest(p):
    if p == 1:
        return False
    if p == 2:
        return True
    if p % 2 == 0:
        return False
    k = 20
    while k:
        a = random.randint(2, p - 2)
        if GCD(a, p) != 1:
            return False
        if PowMod(a, p - 1, p) != 1:
            return False
        k -= 1
    return True

print(FermatTest(int(input())))

Fermat测试的复杂度还能被再次降低。

若$p$为奇数，则有$p-1=2^st$，其中$t$为奇数，则有以下分解式：

$n^{p-1}-1=(n^t-1)(n^t+1)(n^{2t}+1)\times ...\times (n^{s-1}t+1)$

因此，如果有$n^{p-1} \equiv 1(mod p)$，则必有$n^t \equiv 1(mod p)$或者$n^{kt} \equiv -1(mod p)(1 \leqslant k \leqslant s-1)$成立。

故有了Miller-Rabin算法：

• 将$p-1$表示成$2^st$
• 随机取一个数$n(2 \leqslant n \leqslant p-1)$
• 计算$n^t(mod p)$，若其等于$\pm1$，则假定$p$为素数
• 计算$n^{kt}(mod p)(2 \leqslant k \leqslant s-1)$，若其中有一个等于$-1$，则假定$p$为素数。否则$p$为合数
• 重复上述过程$k$次

若测试得到$p$为素数，则准确率为$1-\frac{1}{4^k}$.代码如下：

 

import random


def GCD(a, b):
    while b:
        a, b = b, a % b
    return a


def PowMod(a, n, p):
    r, t = 1, a
    while n:
        if n & 1:
            r = (r * t) % p
        t = (t * t) % p
        n >>= 1
    return r


def MillerRabin(p):
    if p == 1:
        return False
    if p == 2:
        return True
    if p % 2 == 0:
        return False
    s, t = 0, p - 1
    while t % 2 == 0:
        t >>= 1
        s += 1
    k = 10
    while k:
        a = random.randint(2, p - 2)
        if GCD(a, p) != 1:
            return False
        m = PowMod(a, t, p)
        npass = 1
        if m == 1 or m == p - 1:
            npass = 0
        q = s - 1
        while npass and q:
            m = (m * m) % p
            if m == p - 1:
                npass = 0
            q -= 1
        if npass:
            return False
        k -= 1
    return True

print(MillerRabin(int(input())))