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

PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 343 / 分数数列

正の整数 k に対して分数 xi/yi による有限数列 ai は次のように定義される。

a1 = 1/k
ai = (x(i-1)+1)/(y(i-1)-1) [i>1 で約分可能なときは約分する]

ai がある整数 n になったとき数列はそこで終了とし,f(k) = n とする。

たとえば k = 20 のときは次のようになる。

1/20 → 2/19 → 3/18 = 1/6 → 2/5 → 3/4 → 4/3 → 5/2 → 6/1 = 6

f(20) = 6 である。

同様に f(1) = 1, f(2) = 2, f(3) = 1 であり,
Σ[k=1,100] f(k^3) = 118937 である。

Σ[k=1, 2*10^6] f(k^3) を求めよ。

Problem 343 - Project Euler

とりあえず答えを出す

問題文中の例が非常にいいヒントになっています。

この変換で分母と分子の和は一定。1+20=21 からスタートして 3(小さい方から2つ目の約数)で割って 7=1+6 になります。この後の分数はすべて分母と分子が互いに素で,最終的に6になります。

これをもとにコードを書きました。

  1. 小さい方から2つ目の約数で割っていく。大きい方から2つ目の約数になる
  2. 商が Φ(n)=n-1 をみたす数になったら終了(Φはオイラーのファイ関数)

これで約12分かかりました。

1分におさめる

上のコードを手直しして1分におさめることに成功しました。

  • Φ(n)=n-1 は素数条件と同じ。PrimeQ[n] で代用できます
  • i+1 が何度も出てくるので j=i+1 で書き換える
  • 大きい方から2つ目の約数の求め方を変える

これで約6分。半分になりました。

さて,「小さい方から2つ目の約数」とは「最小の素数の約数」です。これでどんどん割っていって最後に残る数は「最大の素数の約数」。実は While で繰り返し割る必要はありません。FactorInteger で素因数分解して,最後の項を拾うだけですみます。

また,k^3+1 をあらかじめ因数分解しておくと大きな数を扱わなくてすみます。このように直したところ,30秒で答えが出ました。