砂場で遊ぼう

c++/python/mathematica などの練習帳

トップ > mathematica > Project Euler 175 / ある数を2のべき乗の和として表せる方法の数の比

Project Euler 175 / ある数を2のべき乗の和として表せる方法の数の比

mathematica project euler

問題概略

整数 n を 2 のべき乗の和で表すことを考える。ただし各数は高々 2 回しか使ってはいけないものとする。この表し方の数を f(n) とする。 f(0)=1 と定義する。

たとえば 10 には 5 通りの異なる表し方があるので f(10)=5 である。

\begin{align*}
10=1+1+8=1+1+4+4=1+1+2+2+4=2+4+4=2+8
\end{align*}

すべての正の有理数 p/q\ (p>0,\, q>0) に対して f(n)/f(n-1)=p/q となる n が少なくとも一つ存在することが示せる。

たとえば f(n)/f(n-1)=13/17 となる最小の n は 241 であり，その 2 進法表示は 11110001 である。
この 2 進数を上の位から順に読んでいくと「1 が 4 つ，0 が 3 つ，1 が 1 つ」となる。
4, 3, 1 を 241 の「短縮された 2 進法表示」と呼ぶ。

f(n)/f(n-1)=123456789/987654321 となる最小の n を短縮された 2 進法表示で求めよ。

https://projecteuler.net/problem=175

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

カルキン・ウィルフ・ツリー

この f(n) は第 169 問でも出てきました。くわしいことはそちらの解説に譲りますが，次の式が成立します。
f(2m+1)=f(m),\, f(2m)=f(m)+f(m-1)

mysandbox.hatenadiary.com

これを使うと f(n)/f(n-1) の引数を減らせます。

\frac{f(2m+1)}{f(2m)}=\frac{f(m)}{f(m)+f(m-1)},\, \frac{f(2m)}{f(2m-1)}=\frac{f(m)+f(m-1)}{f(m-1)}

なんとなく合同式を連想させる式ですが，実は木として図示すると 2019 年の大阪大学の入試でも出題されたカルキン・ウィルフ・ツリー（Calkin – Wilf tree）の矢印を逆向きにしたものになっています。

en.wikipedia.org

f:id:variee:20211228110844p:plain

既約分数の分母と分子は互いに素なので，図 1 のように分母，分子の大きい方から小さい方を引く操作を繰り返せば最後は必ず 1/1=f(1)/f(0) になります。

その経路を逆向きにたどれば f(1)/f(0) から f(n)/f(n-1) にたどりつきます。ゴールの分子の引数が求める n です。

これを実験で確かめてみましょう。 13/17 から 1/1 までの経路は図 3 のようになります。

f:id:variee:20211228110858p:plain

図 3 の経路を逆向きにたどるとき，分子の引数は次のように変化します。

  • 「 \swarrow」のとき m\to 2m+1 なので 2 進法では「1 ビット左にシフトして 1 を足す」
  • 「 \searrow」のとき m\to 2m なので 2 進法では「1 ビット左にシフト」

これらを単に「 +1」「 +0」であらわしたものが図 4 です。問題文に書いてある通り n=241 です。

n=11110001_{(2)}=241_{(10)}

手計算で解ける

123456789/987654321 に対する経路を描いたら手計算で解けてしまいました。

f:id:variee:20211228110914p:plain

mathematica で解く

手計算で解けてしまいましたが，mathematica でも解くことにします。

これまでの議論からほとんどの移動において「移動回数は分母と分子の割り算の商」「移動した後の値は余り」が正しいことがわかっています。

ただ，分母か分子が1になった後は話しが違います。 x/1 からの経路と 1/x からの経路は次のようになります。

f:id:variee:20211228110935p:plain

  • 図 7 の経路を逆にたどるとき，1 に対して +0 を x-1 回行うので 1\underbrace{00\cdots 0}_{\text{$x-1$個}} になる
  • 図 8 の経路を逆にたどるとき，1 に対して +1 を x-1 回行うので \underbrace{11\cdots 1}_{\text{$x$個}} になる

これをふまえたコードを書きました。up が分子で dw が分母です。
これらが 1 でない間は while ループでまとめて処理して，1 になったら個別に処理しています。計算時間は一瞬でした。

    In[]:= AbsoluteTiming[
        solve[p_, q_] := Module[{res = {}, up, dw},
          {up, dw} = {p, q}/GCD[p, q];
          While[up != 1 && dw != 1,
           If[up > dw,
            AppendTo[res, Quotient[up, dw]]; up = Mod[up, dw],
            AppendTo[res, Quotient[dw, up]]; dw = Mod[dw, up]]];
          If[up == 1,
           AppendTo[res, dw];
           Return@Reverse@res];
          If[dw == 1,
           AppendTo[res, {up - 1, 1}];
           Return@Reverse@Flatten@res]];
        ans = solve[123456789, 987654321];]
       
    Out[]= {0.0000738, None}