読者です 読者をやめる 読者になる 読者になる

Problem 277

Project Euler

Project Eulerの割と新しい問題にチャレンジしてみました。
コラッツ列に関する問題です。
なぜいきなり新しい問題を解いてみたかと言いますと、
Project Eulerの正しい遊び方に気付いてしまったからです。


実はプロジェクトオイラーは、
新しい問題が毎週末に公開され、
公開後何分以内に解けるのか!?
というスピードを競うものだったんですね!


いやはや、気付くのが遅すぎました。
とりあえずメールアドレスを登録したら
「今週末の土曜の夕方18:00に新しい問題○○番が投稿されるよ!」
というメールが届くようになったので、
過去問を粛々と解きつつ、
新しい問題にも果敢に挑戦してみようと思います。
最近の問題は結構難しくて頭を捻るものが多いようです。


その一環としてとりあえずわかりやすい277番を解いてみようかと。
この記事を書いたのは今週の日曜日なのですが、
日曜日の夕方に持て余した情熱を発散させる方法としては
なかなか素晴らしいと思います。
(賛同者がいることは期待してませんw)


また、今回から答えの値を明記するのをやめようと思います。
コード書いてるので意味はあまりないのですが、
答えを見たくない人がもしいた場合に迷惑になるといけませんので。


Problem 277
Problem 277(原文)

整数の修正コラッツ列は値 a[1] から始めて次のようにして得られる:
a[n] が 3 で割り切れるならば a[n]+1 = a[n]/3. これを大きな下ステップ "D" と表す.
a[n] を 3 で割った余りが 1 ならば a[n]+1 = (4a[n] + 2)/3. これを大きな上ステップ "U" と表す.
a[n] を 3 で割った余りが 2 ならば a[n]+1 = (2a[n] - 1)/3. これを小さな下ステップ "d" と表す.
数列は a[n] = 1 となれば終了する.

任意の整数が与えられたとき, ステップの列を書き出すことができる.
例えば a[1]=231 なら,
数列 {a[n]}={231,77,51,17,11,7,10,14,9,3,1} は
ステップ "DdDddUUdDD" に対応する.

もちろん, 同じ列 "DdDddUUdDD...." から始まる列は他にもある.
例えば, a[1]=1004064 なら, 
ステップの列は DdDddUUdDDDdUDUUUdDdUUDDDUdDD である.
実際, 1004064 は, 列 DdDddUUdDD から始まる
最小の可能な a[1] > 10 ** 6 である.

列 "UDDDUdddDDUDDddDdDddDDUDDdUUDd" から始まる
最小の a[1] > 10 ** 15 は何か?


まずはコーディングしたものを記してみます。

import time
_S = time.time()

types = {'D':0, 'U':1, 'd':2}

starts = 'UDDDUdddDDUDDddDdDddDDUDDdUUDd'

# (A[i]n + R[i]) / M[i] is always integer. 0 <= i < len(starts)
A, R, M  = [1], [0], [1]
for c in starts:
    idx = types[c]
    
    if idx == 0:
        A.append(A[-1])
        R.append(R[-1])
    elif idx == 1:
        A.append(4 * A[-1])
        R.append(4 * R[-1] + 2 * M[-1])
    elif idx == 2:
        A.append(2 * A[-1])
        R.append(2 * R[-1] - M[-1])
    M.append(M[-1] * 3)
A = A[1:]
R = R[1:]
M = M[1:]

REM = [{'D':0, 'U':1, 'd':2}[starts[0]]]
for i in xrange(1, len(M)):
    a, r, m = A[i], R[i], M[i]
    remainder = REM[-1]
    while remainder < m:
        if (a * remainder + r) % m == 0:
            break
        remainder += M[i - 1]
    remainder %= m
    REM.append(remainder)

limit = '10 ** 15'

N = eval(limit)
n = (N / M[-1] - 1) * M[-1] + remainder
while n <= N:
    n += M[-1]
if n < 10 ** 20:
    print n
else:
    print '%s + %d' % (limit, n - N)

print '%.3f [sec]' % (time.time() - _S)


難しい問題に直面した時ほど、
「これは実験なんだ…」と自分に言い聞かせながら
ついつい徒然なるままにコーディングしてしまいがちなのですが、
コードをきれいにまとめるのもアレなので
とりあえず汚いままコードをさらしてしまいました。


こんな汚いコードを商品として提供したりなどはいたしませんので、
弊社サービスをご利用の方はどうぞご安心ください(笑)。


答えを求めるまでに 0.000 [sec] しかかからず、
また、問題の 10の6乗 というのを 10の1000000乗 という大きな値にしても
わずか 3.313 [sec] で計算できるので、
速度はまぁ申し分ないかなと思います。