Problem 14
（略）（コラッツの問題で）最も長い数列を作る100万より小さい最初の数は？
http://projecteuler.net/index.php?section=problems&id=14
「コラッツの問題」は「角谷の予想」とも言われます。
適当な自然数を選んで、偶数なら2で割り、奇数なら3倍して1を足します。これを繰り返すと最後には必ず1に辿りつく、というのがこの予想です。例えば、7から始めると、
7 -> 22 -> 11 -> 34 -> 17 -> 52 -> 26 -> 13 ->
40 -> 20 -> 10 -> 5 -> 16 -> 8 -> 4 -> 2 -> 1
これが脅威の問題であるのは、この手の問題はまずループになるからです。社交数もそうですし、類似の問題として考えられる、奇数なら3倍して1を引くなら、
5 -> 14 -> 7 -> 20 -> 10 -> 5
とすぐにループになります。
まずは題意どおりに書いてみましょう。
from itertools import imapdef next(n):
if n % 2 == 0:
return n / 2
else:
return n * 3 + 1def length(n):
l = 0
while n != 1:
n = next(n)
l += 1
return lN = 10 ** 6
def max_length(x, y): return x if x[1] >= y[1] else y
tup = reduce(max_length, imap(lambda n: (n, length(n)), xrange(2, N)))
print tup[0]
しかしこれはとても遅いです。Project Eulerの問題は、素直に書くとだいたい遅いです。この対策として2つの常套手段があります。一つはまとめて解くです。エラトステネスのふるいに代表される方法です。ここでは100万までの配列を用意して、一度通ったらそこまで何歩かを記録しておきます。例えば、
2 -> 1
3 -> 10 -> 5 -> 16 -> 8 -> 4 -> 2
6 -> 3
2は1歩、3は7歩、10は6歩、…、6は8歩だと分かります。
from itertools import imapdef next(n):
if n % 2 == 0:
return n / 2
else:
return n * 3 + 1def run():
a[1] = 0
for n in xrange(2, N):
b = [ ]
while n >= N or a[n] == -1:
b.append(n)
n = next(n)
l = a[n]
for m in reversed(b):
l += 1
if m < N:
a[m] = lM = 10 ** 6
N = M
a = [ -1 ] * N
run()def max_length(x, y): return x if x[1] >= y[1] else y
tup = reduce(max_length, imap(lambda n: (n, a[n]), xrange(2, N)))
print tup[0]
このやり方は比較的速いのですが、一つの数に対してすでに記録した数に到達するまで配列に記録しておかなければならないのが嫌なところです。
Project Eulerのもう一つの常套手段は、逆から解くです。ここでは、
3 -> 10 -> 5 -> 16 -> 8 -> 4 -> 2 -> 1
の代わりに1から逆に辿ります。
1 <- 2 <- 4 <- 8 <- 16 <- 32 <- 64 <- 128 <- …
16 <- 5 <- 10 <- …
64 <- 21 <- 42 <- …
偶数で3で割って1余るとき、枝分かれが発生します。そのときにスタックにその数を入れて、あとで2倍ずつします。こうすればスタックが一つ必要なだけで、まんべんなくステップ数を得ることができます。ただし、記録の範囲を超えたところで辿るのを打ち切ってしまうので、記録されない数が出てきます。その分は正順に辿ります。
from itertools import imapdef next(n):
if n % 2 == 0:
return n / 2
else:
return n * 3 + 1def track_back():
a[2] = 1
a[4] = 2
stack = [ 4 ]
while len(stack) != 0:
n = stack.pop()
l = a[n]
r = n % 3
n *= 2
while n < N:
l += 1
a[n] = l
if r != 0 and n % 3 == 1:
m = (n - 1) / 3
a[m] = l + 1
stack.append(m)
n *= 2def length(n):
l = a[n]
if l == 0:
m = n
counter = 0
while m >= N or a[m] == 0:
m = next(m)
counter += 1
l = counter + a[m]
a[n] = l
return lM = 10 ** 6
N = M
a = [ 0 ] * N
track_back()def max_length(x, y): return x if x[1] >= y[1] else y
tup = reduce(max_length, imap(lambda n: (n, length(n)), xrange(2, N)))
print tup[0]
やってみると、少し速くなりました。
