PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 277 / 修正コラッツ列

整数の修正コラッツ列は値 a1 から始めて次のようにして得られる。

an が 3 で割り切れるならば
a(n+1) = an/3
これを大きな下ステップ "D" と表す。

an を 3 で割った余りが 1 ならば
a(n+1) = (4an + 2)/3
これを大きな上ステップ "U" と表す。

an を 3 で割った余りが 2 ならば
a(n+1) = (2an - 1)/3
これを小さな下ステップ "d" と表す。

数列は an = 1 となれば終了する

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

もちろん, 同じ列 "DdDddUUdDD...." から始まる列は他にもある。例えば a1=1004064 なら, ステップの列は
DdDddUUdDDDdUDUUUdDdUUDDDUdDD である。

実際, 1004064は, 列 DdDddUUdDD から始まる最小の可能な a1 > 10^6 である。

列 "UDDDUdddDDUDDddDdDddDDUDDdUUDd" から始まる最小の a1 > 10^15 は何か?

Problem 277 - Project Euler

a1=231

まずは問題文中に与えられている a1=231, DdDddUUdDD を確かめるためのコードを書きます。

3で割った余りに応じた変換をすると同時に myList に D, U, d を追加していって,
aList="DdDddUUdDD" と比較するというストレートなロジックです。

a1=1004064

a1=1004064, DdDddUUdDDDdUDUUUdDdUUDDDUdDD について調べました。

  • そのままやると85秒
  • D からはじまることに注目。a1 は3の倍数です。これだと28秒
  • Dd からはじまることに注目。a1 は9で割って6余る数です。これだと9秒

a1 の候補をいかに絞り込むかがポイントだとわかります。しぼりこみは可能だし,線型に効きます。

本題

与えられた条件は
UDDDUdddDDUDDddDdDddDDUDDdUUDd

手計算で30段階もの変換をさかのぼるのはつらいので,Mathematica にやってもらいます。30回変換した後の数は

 \dfrac{4194304\, a_1-21110037246199}{205891132094649}

これが整数になる条件を考えます。

 4194304\, a_1 \equiv 21110037246199\pmod{3^{30}}

ModularInverse で 4194304 の逆元を求めて両辺に掛けると

 a_1 \equiv 96521732651065\pmod{3^{30}}

この値からスタートして 3^30 刻みで調べると,すぐに答えが出ます。

ベストな手順は?

私はこのように解いたのですが,よく考えてみると 3^30 で割って
96521732651065 余る数は UDDDUdd... の条件をみたします。単純に3^30を足していって,はじめて 10^15 を超える数が答えです。上のコードを走らせる必要はありませんでした。1行の While 文で事足ります。

また,私は 96521732651065 を求めるのに対話的に ModularInverse を打ち込んだりしましたが,実は Reduce で一発でできます。

3^30 を5回足せばいいこともわかるのがすごい。

以上をまとめると,Mathematica を使った場合のベストな手順はこうだと思います。

  1. 合成関数を求める
  2. Reduce で整数解を求める