线性筛求约数个数

写欧拉计划第十二题的时候，求约数个数，暴力实在太慢，顺便写篇博客记录一下高效方法

改编自：　https://blog.csdn.net/ControlBear/article/details/77527115

d(i) 表示 i 的约数个数

num[i] 表示 i 的最小素因子的个数

prim[i] 表示 第 i 个素数

素数

根据基本算数定理，每一个大于等于2的正整数，都可以被分解成
$$N=p_1^{a_1}{p}_2^{a_2}...p_n^{a_n}$$
其中，p为素数。

线性筛就是每一次把最小素因子筛出。

mark = [0] * N
prim = [0] * N
def initial():
    cnt = 0
    for i in range(2, N):
        if not mark[i]:
            prim[cnt] = i
            cnt += 1
        for j in range(0, cnt):
            if i * prim[j] < N:
                mark[i * prim[j]] = 1
                if not i % prim[j]:
                    break

其中，mark，为已知质数的倍数。

##　约数

算数基本定理中，根据拆分后的素因子的指数，我们可以求出每个Ｎ的约数的个数
$$d(N)=(1+a_1)(1+a_2)...(1+a_n)$$
根据这个式子，可以用线性筛去筛去当前Ｎ的约数个数。

筛的过程中，我们需要保存最下素因子的个数。

①当前数是素数

当前$d(N)=(1+1)=2$, 因为素数只有一个素因子（它本身），并且指数为１，

最小素因子个数num[N] = 1

② $i\%prim[j]!=0$

i中不包含 prim[j]这个素因子，然而在i*prim[j]中，包含了一个prim[j]

可以从前面得到i的所有约数个数 $(1+a_1)(1+a_2)...(1+a_n)$

然后补上 prim[j]的个数$(1+a_1)(a+a_2)...(1+a_n)\ast (1+1)$

所以最后得到 $d(i\ast prim(j))=d(i)\ast d(prim[j])$

而且由于当前的 prim[j] 必然是 i * prim[j]的最小素因子，要记录下这个最小素因子个数

所以 num[i * prim[j]] = 1

③i % prim[j] == 0

i中必然包含了至少一个prim[j], 而且prim[j]也必然是i的最小素因子。

而i * prim[j]比起i则是多了一个最小素因子个数，即$1+a_1$

那么 i * prim[j]的约数个数应该是$(1+a_1+1)(1+a_2)...(1+a_n)$

之后，i的最小素因子个数为num[i], 而d(i)中已经包含了$(1+a_1)(1+a_2)...(1+a_n)$

这时可以除去第一项 $1+a_1$，然后乘以$1+a_1+1$，就可得到d(i * prim[j])的约数的个数

$d(i\ast prim[j])=d(i)/(num[i]+1)\ast (num[i]+2)$

当前 num[i * prim[j]] = num[i] +1， 继续保存当前最小素因子个数

N = 20
d = [0] * N
prim = [0] * N
num = [0] * N
mark = [0] * N
def initial():
    cnt = 0
    d[1] = 1
    for i in range(2, N):
        if not mark[i]:
            prim[cnt] = i
            cnt += 1
            num[i] = 1
            d[i] = 2
        for j in range(0, cnt):
            if i * prim[j] < N:
                mark[i * prim[j]] = 1
                if not i % prim[j]:
                    num[i * prim[j]] = num[i] + 1
                    d[i * prim[j]] = d[i] / (num[i] + 1) * (num[i*prim[j]]+ 1)
                    break
                d[i * prim[j]] = d[i] * d[prim[j]]
                num[i * prim[j]] = 1