約数の個数を求めるアルゴリズムと約数が多い数である高度合成数についてのメモ
約数の個数の求め方
実際に割ってみる方法
ある数\(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
