唯物是真 @Scaled_Wurm

プログラミング(主にPython2.7)とか機械学習とか

約数の個数の求め方と高度合成数(Highly composite number)のリスト(10^8まで)

約数の個数を求めるアルゴリズムと約数が多い数である高度合成数についてのメモ

約数の個数の求め方

実際に割ってみる方法

ある数\(n\)の約数の数を求めます
\(\sqrt n\)より大きい約数は\(n\)自身しかないので、\(1\)から\(\sqrt n\)まで割れば十分です
割り切れたときにちょうど平方根でなかったら2つ分カウントすればよいです
(最悪)計算量は\(O(\sqrt n)\)ぐらい(注、この記事にかかれている計算量の理解はゆるふわです)
実際に割っているのでメモリーの量が足りれば約数のリストも作れます

import math
def num_divisors_trial_division(n):
    num = 0
    for i in xrange(1, int(math.sqrt(n) + 1e-9) + 1):
        if n % i == 0:
            num += 1
            if n != i**2:
                num += 1
    
    return num

素因数分解を使う方法

素因数分解が済んでいれば、約数の個数は以下のような式で計算できます
ある数\(n\)が以下のように素因数分解できるとき
$$n = \prod_{i=1}^r {p_i}^{a_i}$$
約数の個数\(d(n)\)は以下の式で求められます
$$d(n) = \prod_{i=1}^r (a_i + 1)$$

素因数の数を持ったdictが与えられた場合以下のようなコードで計算できます

def num_divisors(prime_count):
    num = 1
    
    for prime, count in prime_count.iteritems():
        num *= (count + 1)
    
    return num

前に以下の記事でも書きましたがエラトステネスのふるい風のアルゴリズムを使って素因数分解ができるので、あらかじめ表を作っておけば効率的に計算できます
sucrose.hatenablog.com
ところで素因数分解された形から約数をすべて復元するときの計算量ってどれぐらいなんでしょう
下の記事によると約数の個数は\(\log d(n) \leq O \left (\frac{\log n}{\log \log n} \right )\)らしいので変形すると\(d(n) \leq O(n^{\frac{1}{\log \log n}})\)で、素因数から約数を復元する時について逆に約数が揃っている状態から考えると1種類素因数が減るごとに少なくとも\(\frac{1}{2}\)ずつ減っていくので等比級数分のループが動くと考えると、計算量は同じく\(O(n^{\frac{1}{\log \log n}})\)ぐらいでよいのかな?
integers.hatenablog.com

エラトステネスのふるい風のアルゴリズム

上でエラトステネスのふるい風の素因数分解の話をしましたが、エラトステネスのふるい風のアルゴリズムを使って直接的に約数の数を求めることもできます(これも上の記事に書いてあります)

これを使うと非常に簡単かつ効率的に\(1\)から\(n\)までの数の約数の個数を求めることができます
具体的には、長さ\(n\)の配列を用意して\(1\)の倍数に\(+1\)、\(2\)の倍数に\(+1\)、\(3\)の倍数に\(+1\)、\(4\)の倍数に\(+1\)、……と割り切れる数に足していくだけです

以下のように簡潔なコードになります

def num_divisors_table(n):
    table = [0] * (n + 1)
    
    for i in xrange(1, n + 1):
        for j in xrange(i, n + 1, i):
            table[j] += 1
    
    return table

これの計算量は、各整数\(i\)について\(\frac{n}{i}\)回ずつ処理を行うので調和級数の和の\(n\)倍になるので次のようになります
$$O(\sum_{i=1}^n \frac{n}{i} = \frac{n}{1} + \frac{n}{2} + \frac{n}{3} + \frac{n}{4} + \frac{n}{5} + \dots + \frac{n}{n}) = O(n \log n)$$
これもメモリーが許せば一緒に約数のリストを作ることもできます

速度比較

上に書いた3つの方法の速度を比較します(Pythonなので大きな数に対してやるには遅いです)
\(1\)から\(10^6\)までの数について約数の個数を求めるのにかかった秒数の10回の平均を取りました
結果は以下のようになりました

実際に割る方法 素因数分解による方法 エラトステネスのふるい風の方法
34秒 3.7秒 1.4秒

計算量は実際に割る方法が\(O(n\sqrt n)\)で残りの2つが\(O(n \log n)\)ぐらいのはず

ソースコード
import timeit

trial_division_setup = """
import math
def num_divisors_trial_division(n):
    num = 0
    for i in xrange(1, int(math.sqrt(n) + 1e-9) + 1):
        if n % i == 0:
            num += 1
            if n != i**2:
                num += 1
    
    return num

def main(limit):
    num_of_divisors = []
    for i in xrange(1, limit + 1):
        num_of_divisors.append(num_divisors_trial_division(i))
"""

prime_factorization_setup = """
import collections

def prime_factor_table(n):
    table = [0] * (n + 1)
    
    for i in xrange(2, n + 1):
        if table[i] == 0:
            for j in xrange(i + i, n + 1, i):
                table[j] = i
    
    return table

def prime_factor(n, prime_factor_table):
    prime_count = collections.defaultdict(int)
    
    while prime_factor_table[n] != 0:
        prime_count[prime_factor_table[n]] += 1
        n /= prime_factor_table[n]
    
    if n > 1:
        prime_count[n] += 1
    
    return prime_count

def num_divisors(prime_count):
    num = 1
    
    for prime, count in prime_count.iteritems():
        num *= (count + 1)
    
    return num

def main(limit):
    table = prime_factor_table(limit)

    num_of_divisors = []
    for i in xrange(1, limit + 1):
        factor_count = prime_factor(i, table)
        num_of_divisors.append(num_divisors(prime_factor(i, table)))
"""

eratosthenes_setup = """
def num_divisors_table(n):
    table = [0] * (n + 1)
    
    for i in xrange(1, n + 1):
        for j in xrange(i, n + 1, i):
            table[j] += 1
    
    return table

def main(limit):
    num_divisors_table(limit)
"""


print 'trial_division:', timeit.timeit('main(10**6)', trial_division_setup, number=10)
print 'prime_factorization:', timeit.timeit('main(10**6)', prime_factorization_setup, number=10)
print 'eratosthenes:', timeit.timeit('main(10**6)', eratosthenes_setup, number=10)

高度合成数

1以上の整数を考えた時、自分より小さい数よりも約数の数が大きい物を高度合成数と呼ぶらしいです

高度合成数のリスト

というわけで実際に高度合成数を求めたので以下に書いておきます(ソースコードは一番最後に載せました)

高度合成数は小さい順に素数をかけあわせた形で、素数の指数がどんどん減っていくようになっています

何番目か 高度合成数 約数の個数 素因数分解した形
1 1 1 \(1\)
2 2 2 \(2^1\)
3 4 3 \(2^2\)
4 6 4 \(2^1 \times 3^1\)
5 12 6 \(2^2 \times 3^1\)
6 24 8 \(2^3 \times 3^1\)
7 36 9 \(2^2 \times 3^2\)
8 48 10 \(2^4 \times 3^1\)
9 60 12 \(2^2 \times 3^1 \times 5^1\)
10 120 16 \(2^3 \times 3^1 \times 5^1\)
11 180 18 \(2^2 \times 3^2 \times 5^1\)
12 240 20 \(2^4 \times 3^1 \times 5^1\)
13 360 24 \(2^3 \times 3^2 \times 5^1\)
14 720 30 \(2^4 \times 3^2 \times 5^1\)
15 840 32 \(2^3 \times 3^1 \times 5^1 \times 7^1\)
16 1260 36 \(2^2 \times 3^2 \times 5^1 \times 7^1\)
17 1680 40 \(2^4 \times 3^1 \times 5^1 \times 7^1\)
18 2520 48 \(2^3 \times 3^2 \times 5^1 \times 7^1\)
19 5040 60 \(2^4 \times 3^2 \times 5^1 \times 7^1\)
20 7560 64 \(2^3 \times 3^3 \times 5^1 \times 7^1\)
21 10080 72 \(2^5 \times 3^2 \times 5^1 \times 7^1\)
22 15120 80 \(2^4 \times 3^3 \times 5^1 \times 7^1\)
23 20160 84 \(2^6 \times 3^2 \times 5^1 \times 7^1\)
24 25200 90 \(2^4 \times 3^2 \times 5^2 \times 7^1\)
25 27720 96 \(2^3 \times 3^2 \times 5^1 \times 7^1 \times 11^1\)
26 45360 100 \(2^4 \times 3^4 \times 5^1 \times 7^1\)
27 50400 108 \(2^5 \times 3^2 \times 5^2 \times 7^1\)
28 55440 120 \(2^4 \times 3^2 \times 5^1 \times 7^1 \times 11^1\)
29 83160 128 \(2^3 \times 3^3 \times 5^1 \times 7^1 \times 11^1\)
30 110880 144 \(2^5 \times 3^2 \times 5^1 \times 7^1 \times 11^1\)
31 166320 160 \(2^4 \times 3^3 \times 5^1 \times 7^1 \times 11^1\)
32 221760 168 \(2^6 \times 3^2 \times 5^1 \times 7^1 \times 11^1\)
33 277200 180 \(2^4 \times 3^2 \times 5^2 \times 7^1 \times 11^1\)
34 332640 192 \(2^5 \times 3^3 \times 5^1 \times 7^1 \times 11^1\)
35 498960 200 \(2^4 \times 3^4 \times 5^1 \times 7^1 \times 11^1\)
36 554400 216 \(2^5 \times 3^2 \times 5^2 \times 7^1 \times 11^1\)
37 665280 224 \(2^6 \times 3^3 \times 5^1 \times 7^1 \times 11^1\)
38 720720 240 \(2^4 \times 3^2 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
39 1081080 256 \(2^3 \times 3^3 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
40 1441440 288 \(2^5 \times 3^2 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
41 2162160 320 \(2^4 \times 3^3 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
42 2882880 336 \(2^6 \times 3^2 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
43 3603600 360 \(2^4 \times 3^2 \times 5^2 \times 7^1 \times 11^1 \times 13^1\)
44 4324320 384 \(2^5 \times 3^3 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
45 6486480 400 \(2^4 \times 3^4 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
46 7207200 432 \(2^5 \times 3^2 \times 5^2 \times 7^1 \times 11^1 \times 13^1\)
47 8648640 448 \(2^6 \times 3^3 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
48 10810800 480 \(2^4 \times 3^3 \times 5^2 \times 7^1 \times 11^1 \times 13^1\)
49 14414400 504 \(2^6 \times 3^2 \times 5^2 \times 7^1 \times 11^1 \times 13^1\)
50 17297280 512 \(2^7 \times 3^3 \times 5^1 \times 7^1 \times 11^1 \times 13^1\)
51 21621600 576 \(2^5 \times 3^3 \times 5^2 \times 7^1 \times 11^1 \times 13^1\)
52 32432400 600 \(2^4 \times 3^4 \times 5^2 \times 7^1 \times 11^1 \times 13^1\)
53 36756720 640 \(2^4 \times 3^3 \times 5^1 \times 7^1 \times 11^1 \times 13^1\times 17^1\)
54 43243200 672 \(2^6 \times 3^3 \times 5^2 \times 7^1 \times 11^1 \times 13^1\)
55 61261200 720 \(2^4 \times 3^2 \times 5^2 \times 7^1 \times 11^1 \times 13^1\times 17^1\)
56 73513440 768 \(2^5 \times 3^3 \times 5^1 \times 7^1 \times 11^1 \times 13^1\times 17^1\)
ソースコード

ついでに素因数分解もしたかったので、エラトステネスのふるい風の素因数分解を行って約数の数を求めました(たぶん先に約数の個数を求めてから高度合成数だけ素因数分解したほうが効率的な気がしますが)

import collections

def prime_factor_table(n):
    table = [0] * (n + 1)
    
    for i in xrange(2, n + 1):
        if table[i] == 0:
            for j in xrange(i + i, n + 1, i):
                table[j] = i
    
    return table

def prime_factor(n, prime_factor_table):
    prime_count = collections.defaultdict(int)
    
    while prime_factor_table[n] != 0:
        prime_count[prime_factor_table[n]] += 1
        n /= prime_factor_table[n]
    
    if n > 1:
        prime_count[n] += 1
    
    return prime_count

def num_divisors(prime_count):
    num = 1
    
    for prime, count in prime_count.iteritems():
        num *= (count + 1)
    
    return num

LIMIT = 10**8

table = prime_factor_table(LIMIT)

maximum = 0
count = 0
for i in xrange(1, LIMIT + 1):
    factor_count = prime_factor(i, table)
    num = num_divisors(prime_factor(i, table))
    if num > maximum:
        count += 1
        str_factor = r' \times '.join(['{}^{}'.format(k, factor_count[k]) for k in sorted(factor_count)])
        print r'|{}|{}|{}|\({}\)|'.format(count, i, num, str_factor)
        maximum = num
-->