素因数分解のアルゴリズムは数多くありますが、ここでは2から連続した整数すべての素因数分解を高速に行う方法を紹介します。
ステップ １． 最小素因数のリストを作る
まず、エラトステネスの篩的なアルゴリズムで最小素因数のリストを作ります。最初の素因数が見つかったらそれ以上は何もしないという感じです。
def smallest_factors(N):
sfl = [0]*(N+1)
for p in range(2, N+1):
if sfl[p]: continue # not prime
for pp in range(p, N+1, p):
if sfl[pp]: continue
sfl[pp] = p
return sfl
N = 10
sfcts = smallest_factors(N)
print([(i,sfcts[i]) for i in range(N+1)])
# [(0, 0), (1, 0), (2, 2), (3, 3), (4, 2), (5, 5), (6, 2), (7, 7), (8, 2), (9, 3), (10, 2)]
ステップ ２．nを最小素因数で1になるまで割っていく
最小素因数のリストはべき乗の情報もなく$8$も$10$も$2$となります。こんなものが役に立つのかと思いますが。これを使って2から10までの素因数分解をsympyのfactorintと同様のDICT形式のリストにします。
fctlist = []
for n in range(2, N+1):
fct = dict() # factorint list
n1 = n
while n1>1:
d = sfcts[n1]
pw = fct.get(d,0) # 0 if no entry for d
fct[d] = pw+1
n1 //= d
fctlist.append(fct)
print(fctlist)
# [{2: 1}, {3: 1}, {2: 2}, {5: 1}, {2: 1, 3: 1}, {7: 1}, {2: 3}, {3: 2}, {2: 1, 5: 1}]
ちゃんと求まってそうです。何故そうなるのかと言うと、以下のように$n$を最小素因数で1になるまで次々に割って、その素因数のべき乗を１つ増やして行くと、素因数分解となるという、割と当たり前の仕組みです。
10 / 2 \to 5/5 \to 1 \\
従って10 = 2\times 5
【応用】[2..N]の整数のそれそれの階乗の素因数分解
このアルゴリズムの面白いところは、素因数を追加しながら計算するので階乗の場合には、前のコードで毎回DICT変数を初期化していたのをやめるだけで漸化式的に階乗の素因数分解を求めることが出来ます。
fct = dict() # **** moved outside of the loop
for n in range(2, N+1):
n1 = n
while n1>1:
d = sfcts[n1]
pw = fct.get(d,0) # 0 if no entry for d
fct[d] = pw+1
n1 //= d
print(f"{n}!:", fct)
# 2!: {2: 1}
# 3!: {2: 1, 3: 1}
# 4!: {2: 3, 3: 1}
# 5!: {2: 3, 3: 1, 5: 1}
# 6!: {2: 4, 3: 2, 5: 1}
# 7!: {2: 4, 3: 2, 5: 1, 7: 1}
# 8!: {2: 7, 3: 2, 5: 1, 7: 1}
# 9!: {2: 7, 3: 4, 5: 1, 7: 1}
# 10!: {2: 8, 3: 4, 5: 2, 7: 1}
階乗はnが大きくなると非常に値が大きくなってしまうのですが、このように素因数分解の形だと扱いやすくなります。「大きな整数の階乗の剰余を求める」で紹介しているように剰余計算もpow関数が使えます。
（開発環境：Google Colab）
この考え方はProject Euler、Problem 675: 2**w(n)を解くのに役に立ちます。