Project EulerのProblem 12は三角数を昇順に列挙してその約数の個数が500より大きい最初の数を得ていましたが、Eulerianならきっと約数の個数が500より大きい自然数を昇順に列挙して最初に三角数になる数を見つけようとするでしょう。しかし、これが難しくてなかなか出来ませんでした。今回やっと方法をみつけたので紹介します。
まず、簡単のために約数の個数は500でなく50とします。さて、約数の個数が50を超える最も小さな自然数は、
5040 = 24・32・5・7
です。約数の個数は指数の組合せのみで決まって、
(4 + 1)(2 + 1)(1 + 1)(1 + 1) = 60
ですね。上の素因数分解をよく見てみましょう。2つ特徴があります。1つ目は素数が小さい方から順に使われていて飛びが無いことです。これは当然です。例えば、
24・32・5・11
と上と比べたら上のほうが小さいことは明らかでしょう。2つ目は指数が左から降順であることです。これも、
22・34・5・7
と比べれば当然であることがわかるでしょう。
このように指数の組合せが決まれば最小の数が決まることがわかります。わかりやすいように今後は指数の組合せを次のように書きましょう。
(4, 2, 1, 1)
さて、約数の個数が50より大きい指数の組合せを昇順に列挙するのですが、(1)からはじめて指数を一つずつ加えていきましょう。どのように加えるのがいいでしょう。例えば、
(4, 2, 1)
に1を加えた指数の組合せには、
(5, 2, 1), (4, 3, 1), (4, 2, 2), (4, 2, 1, 1)
があります。どれを選択するのがよいでしょう。うまくやらないと重複して生成してしまいます。ここは逆に考えてみます。すなわち、(5, 2, 1)はどの指数の組合せに1を加えたらいいでしょう。候補は、
(4, 2, 1), (5, 1, 1), (5, 2)
です。当然一番左が最も大きいので、一番左を選択します。ということは、基本的に一番左に指数を一つ加えるのがよいということになります。しかし、例えば、(4, 4, 2)はこの方法では生成することはできません。これは(4, 3, 2)から生成することになります。一番左とその隣が一つ差なら、一番左とその隣に指数を加えることができます。
(4, 2, 1) => (5, 2, 1)
(4, 3, 2) => (5, 3, 2), (4, 4, 2)
(4, 4, 3) => (5, 4, 3), (4, 4, 4)
(1, 1, 1) => (2, 1, 1), (1, 1, 1, 1)
最後の1が続くときのみ項数を増やすことができます。値と指数の組合せをタプルしにて、PriorityQueueに投入すれば、値が小さい順に組合せをキューから取り出せます。取り出したら上の方法でまた組合せをキューに投入すれば、いつか約数の個数が50を超える最小の数が出てきます。やっとファーストステップが終わりです。
指数の組合せが出てきたら、次はその並び替えを考えます。隣同士が左＞右になっていたら入れ替えることができます。(3, 2, 1)なら下のようになります。
(3, 2, 1) => (2, 3, 1), (3, 1, 2)
そして次に、
(2, 3, 1) => (2, 1, 3)
(3, 1, 2) => (1, 3, 2)
このようにしていくと全ての並びを得ることができます。なぜかというと、これはバブルソートの逆だからです。どんな並びも隣同士で左＜右になっていたら交換すればいつか降順に並びます。ただ、このやり方は重複が避けられないので、そこはsetを使います。
さきほどと同じようにPriorityQueueを使って、約数の個数が50より大きい組合せが出てきたら、入れ替えてまたPriorityQueueに投入します。こうして、並びが降順でないものも含めて約数の個数が50より大きい数を昇順で列挙することができます。
最後に、それぞれの並びに対して右にずらします。(1, 3, 2)に対して、(1, 0, 3, 0, 2)です。これは最初のステップと同じように考えられます。まず、右端を一つ右へずらすことができます。
(1, 0, 3, 0, 2) => (1, 0, 3, 0, 0, 2)
それから、右から最初の0の左が0でなければそこも右へずらすことができます。
(1, 0, 3, 0, 2) => (1, 0, 0, 3, 2)
これで約数の個数が50より大きいすべての自然数を昇順に並べることができます。
以前のコードと速度を比較します。
| 前回 | 今回 | |
|---|---|---|
| 500 | 0.08s | 0.1s |
| 5000 | 15.6s | 1.2s |
| 50000 | - | 29.5s |
前回のコードは50000だと数時間かかりそうです。
from itertools import * from Queue import PriorityQueue from math import sqrt def is_prime(n): return all(n % p != 0 for p in takewhile(lambda p: p * p <= n, ps)) ps = [ 2 ] def primes(k): if k < len(ps): return ps[k] for n in count(ps[-1] + 1): if is_prime(n): ps.append(n) if len(ps) > k: break return ps[k] def value(es): return reduce(lambda x, (p, e): x * p ** e, izip(imap(primes, count()), es), 1) def num_divs(es): return reduce(lambda x, e: x * (e + 1), es, 1) def gen_greater_combinations(N): pq = PriorityQueue() pq.put((2, (1,))) while not pq.empty(): v0, es0 = pq.get() if num_divs(es0) > N: yield v0, es0 es1 = (es0[0] + 1,) + es0[1:] pq.put((value(es1), es1)) for k in xrange(1, len(es0)): if es0[k-1] == es0[k] + 1: es2 = es0[:k] + (es0[k] + 1,) + es0[k+1:] pq.put((value(es2), es2)) break elif es0[k-1] > es0[k]: break else: if es0[-1] == 1: es2 = es0 + (1,) pq.put((value(es2), es2)) def is_descending(es): return all(e1 >= e2 for e1, e2 in izip(es, es[1:])) def gen_greater_permutations(N): pq = PriorityQueue() g = gen_greater_combinations(N) set_used = set() pq.put(next(g)) while not pq.empty(): v, es0 = pq.get() yield v, es0 for k in xrange(len(es0) - 1): if es0[k] > es0[k+1]: es1 = es0[:k] + (es0[k+1], es0[k]) + es0[k+2:] v1 = value(es1) if v1 not in set_used: pq.put((v1, es1)) set_used.add(v1) if is_descending(es0): try: pq.put(next(g)) except StopIteration: pass def right_isolated_zero_pos(es): for k in xrange(len(es) - 2, 0, -1): if es[k] == 0: if es[k-1] > 0: return k else: return -1 return -1 def gen_greater_numbers(N): pq = PriorityQueue() g = gen_greater_permutations(N) pq.put(next(g)) while not pq.empty(): v, es0 = pq.get() yield v es1 = es0[:-1] + (0, es0[-1]) pq.put((value(es1), es1)) pos = right_isolated_zero_pos(es0) if pos > 0: es2 = es0[:pos-1] + (0, es0[pos-1]) + es0[pos+1:] pq.put((value(es2), es2)) if all(e > 0 for e in es0): try: pq.put(next(g)) except StopIteration: pass def is_triangle(n): m = int(sqrt(n * 2)) return n == m * (m + 1) / 2 N = 50000 print next(ifilter(is_triangle, gen_greater_numbers(N)))
