Project Eulerでしばしば、○○を満たすN未満の最大の整数を求めよ、といった問題を見かけます。例えば、Problem 26です。この問題を読み替えると、「10がFpの生成元になるN未満の最大の素数を求めよ」とほぼなります。こういう問題は、N-1から順番に整数が題意を満たすかどうか調べて最初に満たすものが答えになるだけなので簡単です。
これは1次元だから順番に並べるのが簡単なのです。では2次元以上の場合はどうすればよいでしょうか。
2次元の場合
この問題はN未満の題意を満たす最大の数を求める問題なので、3桁同士を掛けた数を大きい順に並べればよいだけです。しかし、並べるのは2次元なので難しいです。
999 * 999, 999 * 998, 998 * 998, 999 * 997, ...
たぶんこれくらいはわかると思います。ではその次はどうでしょう。難しいですね。
これを簡単に実現するのがPriorityQueueです。少し遊んでみましょう。
from Queue import PriorityQueue pq = PriorityQueue() pq.put(3) pq.put(1) pq.put(4) print pq.get() # 1 print pq.get() # 3 print pq.get() # 4
このようにPythonのPriorityQueueは投入の順番に関わらずQueueの中で最も小さい要素が取り出されます。
さて、上の問題はどうしたらいいでしょう。わかりやすいように3桁ではなく1桁で考えます。まず、最大は9 * 9なのでこれを投入します。小さい順に取り出されるのでマイナスをつけます。そして、積の両項もほしいので、(-9 * 9, 9, 9)とします。一般にx * yのとき(-x * y, x, y)を投入します。
pq = PriorityQueue() pq.put((-9 * 9, 9, 9))
問題は、これを取り出した後です。PriorityQueueは「次」が大事なのです。
p, x, y = pq.get()
pに-9 * 9、xに9、yに9が代入されます。次は隣へ行けるとします。隣とはここではxかyのどちらかが同じでどちらかが1違うxとyのペアとします。それは、9 * 8です。8 * 9は9 * 8と区別をつける必要がないので要りません。常にx >= yとします。
pq.put((-9 * 8, 9, 8))
次に取り出されるのは9 * 8です。この次は8 * 8と9 * 7が隣で両方とも投入します。
次に取り出されるのは8 * 8です。この次は8 * 7です。しかし、9 * 7の次にも8 * 7があります。ここがPriorityQueueを使うときの難しさの本質です。
9 * 9 ↓ 9 * 8 → 8 * 8 ↓ ↓ 9 * 7 → 8 * 7
重複があると効率が悪くなったりする可能性があるので、それをなくすためにルールを設けます。それはxが9のときだけ下にいけるというものです。そうすると次のような図になります。
9 * 9 ↓ 9 * 8 → 8 * 8 ↓ 9 * 7 → 8 * 7 → 7 * 7 ↓ 9 * 6 → 8 * 6 → 7 * 6 → 6 * 6 ↓
このように重複無しに全ての要素を辿ることができるようにすることがPriorityQueueを使うということの本質です。
これを元に3桁の自然数の積を降順に並べるコードを書くと次のようになります。
from itertools import islice from Queue import PriorityQueue # このジェネレータは汎用的に使える def gen_sorted(init, nexts): pq = PriorityQueue() pq.put(init) while not pq.empty(): z0 = pq.get() yield z0 for z in nexts(z0): pq.put(z) def gen_products_descendingly(N): def nexts((p, x, y)): # 本当は2桁にならないようにしなければならない if x > y: yield (-(x - 1) * y, x - 1, y) if x == 999: yield (-x * (y - 1), x, y - 1) M = N - 1 return gen_sorted((-M * M, M, M), nexts) N = 1000 print list(islice(((x, y) for p, x, y in gen_products_descendingly(N)), 6))
[(999, 999), (999, 998), (998, 998), (999, 997), (998, 997), (999, 996)]
高度合成数
ラマヌジャンの高度合成数（highly composite number）という数があります。これはたくさんの素数から合成されているので約数の個数が多い数ということです。具体的には、その数より小さいどの自然数よりも約数の個数が多い自然数のことを言います。最初のほうを列挙すると（カッコ内は約数の個数）、
1(1), 2(2), 4(3), 6(4), 12(6), 24(8), 36(9), ...
この数にはいくつかの特徴があります。その前に自然数を図形で表すことを考えます。例えば、60を素因数分解すると、
60 = 22 * 3 * 5
ですが、これを
* ***
と表すことにします。左から2、3、5の個数を表しています。次に、84を素因数分解すると、
84 = 22 * 3 * 7
です。これを図形で表すと、
* ** *
となります。ところで、約数の個数は、60も84もどちらも指数が2, 1, 1だから、
(2 + 1) * (1 + 1) * (1 + 1) = 12
となります。約数の個数は指数で決まるのだから素数はなるべく小さいほうがよいです。つまり、図にしたときに間が空いているような場合は高度合成数にはなりえません。これで高度合成数の候補は指数の組だけで表すことができることになります。60なら
(2, 1, 1)
と表せます。また、間が空いていない場合でも、(1, 2, 1)より(2, 1, 1)のほうが小さいので、(1, 2, 1)は高度合成数になりえません。一般に高度合成数は、
(e1, e2, ..., em)
e1 ≥ e2 ≥ ... ≥ em
の形をしています。
さて、上の形の自然数を高度合成数の候補として、これを小さいものから並べて、今までで一番約数の個数が大きくなっていればそれは高度合成数です。この列挙をPriorityQueueを用いて行います。PriorityQueueに突っ込む要素は、
(自然数, 指数の組, 約数の個数)
とします。初期値は(2, (1,), 2)とします。しかし、この問題は「次」が難しいです。何しろ1次元でも2次元もなく無限次元だからです。例えば、指数を一つ加えて(3, 2, 1)を作りたいとしましょう。(2, 2, 1)でも(3, 1, 1)でも(3, 2)でも(3, 2, 1)になれます。3つのうちからどれを選んだらよいでしょう。
こういうときは、最も効率のよい選択法を選びます。すなわち、自然数がなるべく大きいものから遷移すればPriorityQueueの滞留時間が短くなって効率がよくなります。つまり、なるべく小さな素数で割って得られる自然数から遷移するのがよいということになります。(3, 2, 1)なら
(2, 2, 1) → (3, 2, 1)
となります。しかし、(2, 2, 2, 1)なら2で割ると(1, 2, 2, 1)となって高度合成数の候補になりえないので、
(2, 2, 1, 1) → (2, 2, 2, 1)
とします。これを逆に考えて、次のようなルールにするとうまくいきます。
- 最も左はいつもインクリメントできる
- その指数より左の指数が全て同じでその指数より1大きければインクリメントできる
- 指数が1並びの場合のみ、右に指数1を追加できる
2番目は、(3, 3, 3, 2, 1) → (3, 3, 3, 3, 1)という感じです。
最初のほうだけ図にすると、
(1,) ━ (2,) ━ (3,) ┃ (1, 1) ━ (2, 1) ┳ (3, 2) ┃ ┗ (3, 1) (1, 1, 1)
となります。
以上をコードにして、最初の100個の高度合成数を生成してみました。
from itertools import * from Queue import PriorityQueue def gen_sorted(init, nexts): pq = PriorityQueue() pq.put(init) while not pq.empty(): z0 = pq.get() yield z0 for z in nexts(z0): pq.put(z) def gen_highly_synthized(): ps = [ 2 ] def is_prime(n): return all(n % p != 0 for p in takewhile(lambda p: p * p <= n, ps)) def primes(k): while k >= len(ps): n0 = ps[-1] + 1 p = next(n for n in count(n0) if is_prime(n)) ps.append(p) return ps[k] def nexts((n, es, nd)): e0 = es[0] yield (n * 2, (e0 + 1,) + es[1:], nd / (e0 + 1) * (e0 + 2)) for k in xrange(1, len(es)): if es[k] == e0 - 1: yield (n * primes(k), es[:k] + (e0,) + es[k+1:], nd / e0 * (e0 + 1)) break if es[k] < e0 - 1: break else: if e0 == 1: yield (n * primes(len(es)), es + (1,), nd * 2) max_nd = 1 for n, es, num_divs in gen_sorted((2, (1,), 2), nexts): if num_divs > max_nd: yield n, num_divs max_nd = num_divs print list(islice(gen_highly_synthized(), 100))
