今回は「ハミングの問題 (Hamming's Problem)」を解いてみましょう。
7 以上の素数で割り切れない正の整数は、素因子が 2, 3, 5 しかない自然数のことです。これを「ハミング数 (Hamming Numbers)」といいます。ハミング数は素因数分解したとき、\(2^i \times 3^j \times 5^k \ (i, j, k \geq 0)\) の形式になります。たとえば、100 以下のハミング数は次のようになります。
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100
同様に、素因子が 2, 3, 5, 7 しかない自然数のことを「ハンブル数 (humble number)」といいます。今回は問題を少し変えて、正の整数 n 以下のハミング数とハンブル数をすべて求めるプログラムを Python3 (PyPy3) で作ってみましょう。
一番簡単な方法は、1 から n までの整数列を生成して、そこからハミング数 (またはハンブル数) を取り出していくことです。これを プログラムすると次のようになります。
リスト : ハミングの問題
import time
def check_hamming(n):
while n % 2 == 0: n //= 2
while n % 3 == 0: n //= 3
while n % 5 == 0: n //= 5
return n == 1
def hamming(n):
return [x for x in range(1, n + 1) if check_hamming(x)]
def check_humble(n):
while n % 2 == 0: n //= 2
while n % 3 == 0: n //= 3
while n % 5 == 0: n //= 5
while n % 7 == 0: n //= 7
return n == 1
def humble(n):
return [x for x in range(1, n + 1) if check_humble(x)]
def test0(func, s, e):
for x in range(s, e):
s = time.clock()
print(10 ** x, len(func(10 ** x)))
print(time.clock() - s)
関数 check_hamming(n) は n がハミング数かチェックします。これは 2, 3, 5 だけで割り切れるか試しているだけです。同様に、関数 check_humble(n) は n がハンブル数かチェックします。実行結果は次のようになります。
>>>> hamming(100) [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100] >>>> humble(100) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36, 40, 42, 45, 48, 49, 50, 54, 56, 60, 63, 64, 70, 72, 75, 80, 81, 84, 90, 96, 98, 100] >>>> test0(hamming, 2, 9) 100 34 0.0 1000 86 0.0 10000 175 0.03125 100000 313 0.015625 1000000 507 0.046875 10000000 768 0.453125 100000000 1105 4.453125 >>>> test0(humble, 2, 9) 100 46 0.0 1000 141 0.015625 10000 338 0.03125 100000 694 0.03125 1000000 1273 0.078125 10000000 2155 0.578125 100000000 3427 5.75 実行環境 : PyPy3 (ver 7.3.1), Ubunts 20.04 (WSL1), Intel Core i5-6200U 2.30GHz
プログラムはとても簡単ですが、引数 n の値が大きくなると時間がかかるようになります。n に比べてハミング数の個数は少ないようなので、式 \(2^i \times 3^j \times 5^k \ (i, j, k \geq 0)\) を使ってハミング数を生成したほうがよさそうです。引数 \(n\) に対して \(i, j, k\) の上限値は \(\log_2 n, \log_3 n, \log_5 n\) で求めることができます。たとえば、100000000 の場合は次のようになります。
i : 0 - 26 j : 0 - 16 k : 0 - 11
全体で 27 * 17 * 12 = 5508 個しかありません。この中から 100000000 以下の数を選べばいいわけです。プログラムは次のようになります。
リスト : ハミングの問題 (2)
import math
def hamming2(n):
xs2 = [2 ** x for x in xrange(0, int(math.log(n, 2)) + 1)]
xs3 = [3 ** x for x in xrange(0, int(math.log(n, 3)) + 1)]
xs5 = [5 ** x for x in xrange(0, int(math.log(n, 5)) + 1)]
return sorted([x * y * z for x in xs2 for y in xs3 for z in xs5 if x * y * z <= n])
def humble2(n):
xs2 = [2 ** x for x in range(0, int(math.log(n, 2)) + 1)]
xs3 = [3 ** x for x in range(0, int(math.log(n, 3)) + 1)]
xs5 = [5 ** x for x in range(0, int(math.log(n, 5)) + 1)]
xs7 = [7 ** x for x in range(0, int(math.log(n, 7)) + 1)]
return sorted([a * b * c * d for a in xs2 for b in xs3 for c in xs5 for d in xs7 if a * b * c * d <= n])
2, 3, 5, (7) のべき乗の集合を生成し、その要素を内包表記で掛け合わせて、条件を満たす数値を選択していくだけです。実行結果は次のようになりました。
>>>> test0(hamming2, 8, 12) 100000000 1105 0.03125 1000000000 1530 0.0 10000000000 2053 0.0 100000000000 2683 0.0 >>>> test0(humble2, 8, 12) 100000000 3427 0.015625 1000000000 5194 0.0 10000000000 7575 0.015625 100000000000 10688 0.0
とても速くなりましたね。このほかにも 参考文献 のプログラムが高速で、Python でジェネレータにすることも簡単です。プログラムは次のようになります。
リスト : ハミングの問題 (3)
def hamming3():
hs = []
j2 = j3 = j5 = 0
m2 = m3 = m5 = 1
while True:
m = min(m2, m3, m5)
hs.append(m)
yield m
while m2 <= m:
m2 = 2 * hs[j2]
j2 += 1
while m3 <= m:
m3 = 3 * hs[j3]
j3 += 1
while m5 <= m:
m5 = 5 * hs[j5]
j5 += 1
def humble3():
hs = []
j2 = j3 = j5 = j7 = 0
m2 = m3 = m5 = m7 = 1
while True:
m = min(m2, m3, m5, m7)
hs.append(m)
yield m
while m2 <= m:
m2 = 2 * hs[j2]
j2 += 1
while m3 <= m:
m3 = 3 * hs[j3]
j3 += 1
while m5 <= m:
m5 = 5 * hs[j5]
j5 += 1
while m7 <= m:
m7 = 7 * hs[j7]
j7 += 1
>>>> for i,x in enumerate(hamming3()): .... if i > 100: break .... print(x, end=' ') .... 1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 40 45 48 50 54 60 64 72 75 80 81 90 96 100 108 120 125 128 135 144 150 160 162 180 192 200 216 225 240 243 250 256 270 288 300 320 324 360 375 384 400 405 432 450 480 486 500 512 540 576 600 625 640 648 675 720 729 750 768 800 810 864 900 960 972 1000 1024 1080 1125 1152 1200 1215 1250 1280 1296 1350 1440 1458 1500 1536 1600 >>>> for i,x in enumerate(humble3()): .... if i > 100: break .... print(x, end=' ') .... 1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 20 21 24 25 27 28 30 32 35 36 40 42 45 48 49 50 54 56 60 63 64 70 72 75 80 81 84 90 96 98 100 105 108 112 120 125 126 128 135 140 144 147 150 160 162 168 175 180 189 192 196 200 210 216 224 225 240 243 245 250 252 256 270 280 288 294 300 315 320 324 336 343 350 360 375 378 384 392 400 405 420 432 441 448 450 480
なお、ハミングの問題は遅延ストリーム (遅延リスト) を使って解くことができます。興味のある方は拙作のページをお読みくださいませ。
最後に、Project Euler の Problem 204 (日本語訳) に挑戦してみましょう。この問題は、N 以下の素因数しか持たない正整数を Type N の一般化ハミング数と定義します。すると、ハミング数は Type 5、ハンブル数は Type 7 の一般化ハミング数になります。問題は 109 以下の Type 100 の一般化ハミング数の個数を求めるのですが、本ページでは少し変更して 1012 以下の Type 50 の一般化ハミング数の個数を求めることにします。
プログラムは次のようになります。
リスト : 一般化ハミング数
# 50 以下の素数
ps = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
def general_hamming(n, i, xs, limit):
if len(xs) == i: return 1
m = xs[i]
j = 0
c = 0
while True:
k = n * (m ** j)
if k > limit: break
c += general_hamming(k, i + 1, xs, limit)
j += 1
return c
def solver(xs, limit):
return general_hamming(1, 0, xs, limit)
実際の処理は関数 general_hamming() で行います。引数 n が生成中のハミング数、xs が素因子を格納したリスト、i が乗算する素因子の添字、limit がハミング数の上限値です。最初に i が len(xs) と等しいかチェックします。そうであれば、素因子をすべて乗算して、ハミング数がひとつ完成しました。return で 1 を返します。
次に、乗算する素因子 xs[i] を取り出して変数 m にセットします。while ループの中で n * mj を計算し、変数 k にセットします。k が limit を超えていなければ、次の素因子と乗算するため general_hamming() を再帰呼び出しして、返り値を変数 c に加算します。それから、j を +1 して素因子 m の計算を続けます。最後に return で c を返します。
実行結果を示します。
>>>> solver([2,3,5], 10 ** 11) 2683 >>>> solver([2,3,5], 10 ** 12) 3429 >>>> solver([2,3,5,7], 10 ** 11) 10688 >>>> solver([2,3,5,7], 10 ** 12) 14672 >>>> s = time.time(); print(solver(ps, 10 ** 12)); time.time() - s 9180401 1.9752330780029297
1012 以下のハミング数とハンブル数は高速に求めることができますが、Type 50 の一般化ハミング数になると、ちょっと時間がかかるようになります。興味のある方はいろいろ試してみてください。