素数筛和质因数分解

这里总结一下计算$[1, n]$区间素数的方法。最直接的方式是对区间中的每个数进行枚举，时间复杂度为$O(n\sqrt{n})$，在$n$较大的时候并不好用。这里介绍一下素数筛法，筛法的思路也很直观：对于不超过$n$的每个非负数$p$，删除$2p,3p,4p...$，当处理完所有的数后，最后剩下的就是素数。下面我们从最朴素的筛法开始，再介绍筛法的改进版本。

1. 朴素筛法

将上面筛法的思路直接翻译成代码就是朴素的素数筛法，代码如下...

# 计算[1, n]区间的素数
def prime(n: int) -> List[int]:
    vis = [0] * (n + 1)
    for i in range(2, n + 1):
        for j in range(i * 2, n + 1, i):  # 对2i、3i、4i...进行更新
            vis[j] = 1
    return [i for i in range(1, n + 1) if vis[i] == 0]

内存循环的次数为$\left \lfloor \frac{n}{i} \right \rfloor - 1 < \frac{n}{i}$，循环的总次数小于$\sum_{i=2}^{n}\frac{n}{i}=O(nlgn)$。这个结论来源于欧拉在1734年得到的结果：$\sum_{i=1}^{n}\frac{1}{i}=ln(n+1)+\gamma $，其中欧拉常数$\gamma \approx 0.577218$。

2. 改进筛法

下面来改进上面的朴素筛算法。如果当前的数是合数，就没必要用它作为筛孔了，因为它的因子已经将以它为因子的合数全部筛掉了；另外，内层循环也没必要从$2i$开始，直接从$i^2$开始即可，因为$2i,3i....,(i-1)i$都被$i$之前的素数消除掉了。基于以上两点，改进的代码如下....

# 计算[1, n]区间的素数
def prime(n: int) -> List[int]:
    vis = [0] * (n + 1)
    for i in range(2, n + 1):
        if not vis[i]:  # 在朴素法基础上加了判断
            for j in range(i ** 2, n + 1, i):  # 朴素法为i*2，这里是i**2，注意区别
                vis[j] = 1
    return [i for i in range(1, n + 1) if vis[i] == 0]

这份代码还不是最优的，例如12会被2、3都删除一次，有没有什么办法可以保证每个合数只被删除一次呢？

3. 快速线性筛法

• 唯一性定理：每一个比1大的自然数$N$只能有一种方式分解成质数的乘积。
• 唯一性定理推论：若一个质数$p$是乘积$ab$的因子，则$p$不是$a$的因子就是$b$的因子。（这个只是拓展一下，这里面没用到这个推论）

这里先给出素数分解唯一性定理及其推论，相关的证明《几何原本》或者维基百科上都有，可以自行查阅。

不论是素数还是合数，都可以分解为除1以外的素因子的乘积（素数的话就是其本身），例如$12=2*2*3$。在改进版本的素数筛中，12会被2和3各消除一次，有什么办法可以让12只被消除一次？如果可以让每个合数只被它最小的素因子（除1）消除，就可以保证只被消除一次了，下面先给出代码，然后证明：

1. 当前扫描到的数是合数，这个数肯定已经被消除
2. 每个合数只会被消除一次

# 计算[1, n]区间的素数
def quickPrime(n: int) -> List[int]:
    vis = [0] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if not vis[i]:
            primes.append(i)
        # 这里不论i是素数还是合数都会执行
        for prime in primes:
            if i * prime > n:
                break
            vis[i * prime] = 1
            # 这里是关键
            if i % prime == 0:
                break
    return primes

我们令合数$i=p_1^{r_1}p_2^{r_2}...p_m^{r_m}$，其中$p_1,p_2...p_m$都是素数，且有$p_j < p_k,1\leq j< k\leq m$，假设当前扫描到$i$，$i$肯定已经被$p_1$和$p_1^{r_1-1}p_2^{r_2}...p_m^{r_m}$消掉。不失一般性，我们假设$i$会被$p_2$和$p_1^{r_1}p_2^{r_2-1}...p_m^{r_m}$消掉，从程序中可以看出，$p_1^{r_1}p_2^{r_2-1}...p_m^{r_m}$在执行到$p_1$的时候就会退出，不会再执行到$p_2$，这就保证了每个合数只被消除一次。

4. 质因数分解

改进上面的快速素数筛，在vis数组（对应下面的factors数组）中存的不是一个数是否被访问过，而是存这个数的最小的质因子，就可以递归进行质因子分解。这里有一个预处理过程，如果待分解的数的最大值为10^6，就需要对10^6规模的数据做一个预处理，时间复杂度为O(n)，每个数分解素因子的时间复杂度为O(这个数质因子的个数)。代码如下：

def get_factors(n):  # 返回[2, n]的最小质因子
    factors, primes = [0] * (n + 1), []
    for i in range(2, n + 1):
        if not factors[i]:
            factors[i] = i  # 质数的最小的素因子就是其本身
            primes.append(i)
        for prime in primes:
            if i * prime > n:
                break
            factors[i * prime] = prime  # 每个vis[i]保存的是最小的质因子
            if i % prime == 0:
                break
    return factors


def decomposition(num, factors):  # 注意要满足 len(factors) >= num and num > 1
    ans = []  # 用于保存素因子
    while num > 1:
        temp = factors[num]  # 获取最小的素因子
        ans.append(temp)
        num //= temp
    return ans  # 结果是按照素因子排序的

还有一个分解质因子的方法，假设待分解的数的最大值为10^6，那可以肯定的是这个数一定不会有超过两个大于1000的素因子，用1000以内的素因子分别去除这个数就行了，除完以后判断最后剩下的数是多少，如果大于1，最后剩的这个数就是大于1000的唯一素因子。预处理时间复杂度为O(sqrt(n)), 每个数分解质因子的最坏时间复杂度为O(sqrt(n)以内质数的个数+这个数质因子的个数)，质数的个数可以用质数定理得到。代码如下：

def quickPrime(n):  # 这里用的是快速素数筛法
    vis, primes = [0] * (n + 1), []
    for i in range(2, n + 1):
        if not vis[i]:
            primes.append(i)
        for prime in primes:
            if i * prime > n:
                break
            vis[i * prime] = 1
            if i % prime == 0:
                break
    return primes


def decomposition(num, primes):
    ans = []
    for prime in primes:
        if not num % prime:  # prime是num的素因子
            while not num % prime:  # 可能不只一个prime因子
                num //= prime
                ans.append(prime)
    if num != 1:  # 例如num=10**6，primes是1000以内的质数，最后剩下的num可能是一个大于1000的质数
        ans.append(num)
    return ans  # 结果是按照素因子排序的