线性筛求约数个数

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

改编自：　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 N = p_1^{a_1}p_2^{a_2}...p_n^{a_n} N=p1a1​​p2a2​​...pnan​​
其中，p为素数。

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

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

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

##　约数

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

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

①当前数是素数

当前$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 ) (1+a_1)(1+a_2)...(1+a_n) (1+a1​)(1+a2​)...(1+an​)

然后补上 prim[j]的个数 ( 1 + a 1 ) ( a + a 2 ) . . . ( 1 + a n ) ∗ ( 1 + 1 ) (1+a_1)(a+a_2)...(1+a_n)*(1+1) (1+a1​)(a+a2​)...(1+an​)∗(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 1+a_1 1+a1​

那么 i * prim[j]的约数个数应该是 ( 1 + a 1 + 1 ) ( 1 + a 2 ) . . . ( 1 + a n ) (1+a_1+1)(1+a_2)...(1+a_n) (1+a1​+1)(1+a2​)...(1+an​)

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

这时可以除去第一项 1 + a 1 1+a_1 1+a1​，然后乘以 1 + a 1 + 1 1+a_1+1 1+a1​+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 = 0d[1] = 1for i in range(2, N):if not mark[i]:prim[cnt] = icnt += 1num[i] = 1d[i] = 2for j in range(0, cnt):if i * prim[j] < N:mark[i * prim[j]] = 1if not i % prim[j]:num[i * prim[j]] = num[i] + 1d[i * prim[j]] = d[i] / (num[i] + 1) * (num[i*prim[j]]+ 1)breakd[i * prim[j]] = d[i] * d[prim[j]]num[i * prim[j]] = 1

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！