砂場で遊ぼう

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

Project Euler 143 / 三角形のトリチェリ点

mathematica project euler グラフ ピタゴラス数

問題概略

ABC をすべての内角が 120 度未満の三角形とする。
X を三角形の内点として \mathrm{XA}=p, \mathrm{XC}=q, \mathrm{XB}=r とおく。

フェルマーは「 p+q+r を最小にする X を探す方法はあるか？」とトリチェリに問題を出した。

正三角形 AOB, BNC, AMC が ABC の各辺に構成できるならば，AOB, BNC, AMC に外接する 3 つの円が三角形の内部の 1 点 T で交わることをトリチェリは示した。

さらに，フェルマー点（トリチェリ点，等角中心）と呼ばれる点 T が p+q+r を最小にすることも示した。
p+q+r が最小のとき \mathrm{AN}=\mathrm{BM}=\mathrm{CO}=p+q+r であり，AN, BM, CO が点 T を通ることも示せる。

f:id:variee:20220106203848p:plain

p+q+r が最小で a, b, c, p, q, r がすべて自然数のとき，三角形 ABC をトリチェリ三角形と呼ぶ。
たとえば a=399, b=455, c=511 は p+q+r=784 のトリチェリ三角形である。

トリチェリ三角形について，異なる値をとる p+q+r\ (\leqq 120{,}000) の総和を求めよ。

https://projecteuler.net/problem=143

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

drive.google.com

余弦定理→整数のパラメータ表示

与えられた図において点 T のまわりの角はすべて 120 度です。
\triangle\mathrm{TAB}, \triangle\mathrm{TBC}, \triangle\mathrm{TCA} に余弦定理を使います。

q^2+qr+r^2=a^2,\, p^2+pq+q^2=b^2,\, r^2+rp+p^2=c^2

これをみたす自然数はピタゴラス数の一般項 \left(k(m^2-n^2),\, 2kmn,\, k(m^2+n^2)\right) と同じ要領でパラメータ表示することができます。

q^2+qr+r^2=a^2 を a^2 で割って q/a=x, r/a=y とおきます。 x, y は有理数です。

x^2+xy+y^2=1

この楕円上の格子点として (\pm 1,\, 0) と (0,\, \pm 1) がすぐみつかります。

(-1,\, 0) を通る直線として y=t(x+1) を考えて楕円との交点を求めると， (-1,\, 0) でない方の座標は次のようになります。

(x,\, y)=\left(\frac{1-t^2}{1+t+t^2},\, \frac{2t+t^2}{1+t+t^2}\right)

t=n/m（ m, n は互いに素な自然数）とおきます。

(x,\, y)=\left(\frac{q}{a},\, \frac{r}{a}\right) =\left(\frac{m^2-n^2}{m^2+mn+n^2},\, \frac{2mn+n^2}{m^2+mn+n^2}\right)

(q,\, r,\, a)=(m^2-n^2,\, 2mn+n^2,\, m^2+mn+n^2) がみつかりました。

次はこれらが互いに素な自然数になる条件を求めます。

「 m と n は互いに素」と q>0 からわかる「 m>n」の他に条件がないか mathematica で調べてみました。

  In[]:= ClearAll["Global`*"];
    q[m_, n_] := m^2 - n^2;
    r[m_, n_] := 2 m*n + n^2;
    isBad[{m_, n_}] := 
      m > n && CoprimeQ[m, n] && ! CoprimeQ[q[m, n], r[m, n]];
    Select[Tuples[Range@10, 2], isBad]
  
  Out[]= {{4, 1}, {5, 2}, {7, 1}, {7, 4}, {8, 5}, {10, 1}, {10, 7}}

m-n が 3 の倍数だとまずいみたいですね。これを証明します。

q, r の最大公約数を \langle q,\, r \rangle のようにあらわします。

\langle q,\, r\rangle=\langle m^2-n^2,\, 2mn+n^2\rangle =\langle m^2+2mn,\, 2mn+n^2\rangle=\langle m(m+2n),\, n(2m+n)\rangle

m と n は互いに素なので「 m と n, 2m+n は互いに素」と「 n と m, m+2n は互いに素」がいえます。
上式の最右辺の値は \langle m+2n,\, 2m+n \rangle と同じです。

\langle q,\, r\rangle=\langle m+2n,\, 2m+n \rangle

M=m+2n, N=2m+n とおくと 3m=-M+2N, 3n=2M-N がいえます。
M と N の最大公約数は 3m と 3n の公約数です。1 か 3 ですね。

この 3 を否定すればいいので， q と r が互いに素な条件は M=m+2n, N=2m+n の少なくとも一方が 3 で割りきれないことです。
casework により「 m\not \equiv n\pmod{3}」がわかります。

まとめましょう。

  • q^2+qr+r^2 が平方数になるのは q=k(m^2-n^2), r=k(2mn+n^2) のとき
  • m と n は互いに素な自然数で m>n かつ m\not \equiv n\pmod{3} をみたす
  • k は自然数

整数の組を数える

FindCycle で長さ 3 の閉路を探す


上の条件をみたす 2 つの数を辺で結んでグラフを作ります。

「 p と q」「 q と r」「 r と p」が平方数を与えるとき， p, q, r は長さ 3 の閉路を作ります。

mathematica の組み込み関数 FindCycle でこれを探して， p+q+r\leqq 12\cdot 10^4 をみたすものの和をとると解けるはず……と思いましたが実際には時間がかかりすぎてダメですね。
「 12\cdot 10^3 以下」は 15 秒程度で答えが出ますが，「 12\cdot 10^4 以下」だと固まってしまいます。

  In[]:= AbsoluteTiming[
    lim = 12*10^3;
    g = {};
    For[m = 2, m^2 + 2 m < lim, m++,
     For[n = 1, n < m, n++,
      If[! CoprimeQ[m, n] || Divisible[m - n, 3], Continue];
      p = m^2 - n^2;
      q = 2*m*n + n^2;
      Do[AppendTo[g, k*p <-> k*q], {k, 1, Quotient[lim, p + q]}]]];
    ans = Total@Select[DeleteDuplicates[
        Total /@ VertexList /@ FindCycle[g, {3}, All]], # < lim &]]
   
  Out[]= {15.2798, 251752}

遅い原因は FindCycle と AppendTo です。

  • FindCycle ではデフォルトで閉路を 1 つ探すようになっている。 閉路を全部探すのは現実的ではない
  • リストに新しい要素を追加する AppendTo は遅い

stackoverflow.com

隣接リストを使う

グラフを無向辺の集合として作るのをやめて，隣接リストで管理します。

隣接リストは小さい方から大きい方に辺を張って作ります。
p, q, r (p< q< r) がループを作っていることの確認は「 p\to q, p\to r」→「 q\to r」の 2 段階に分けて行います。

  1. p の隣接リストから大きさ 2 の部分集合 x,\, y\}\ (x< y) を作る。（辺 p\to x, p\to y がある）
  2. x の隣接リストに y が含まれていれば OK（辺 x\to y がある）

他にもいくつか工夫しました。

ループさせない

mathematica はループ処理が苦手です。 (m,\, n) のリストを Table で一気に作りました。
要素数は約 15,000 です。この変更で計算時間はほとんど変わりませんでした。

  lim = 12*10^4;
  mnTable = Flatten[Table[{m, n}, {m, 2, Sqrt[1 + lim] - 1},
     {n, 1, Min[m, Quotient[lim - m^2, 2 m]]}], 1];
  mnTable = 
   Select[mnTable, 
    CoprimeQ[#[[1]], #[[2]]] && ! Divisible[#[[1]] - #[[2]], 3] &];

m, n の範囲は p+q=(m^2-n^2)+(2mn+n^2)=m^2+2mn が \mathrm{lim}=12{,}000 を下回るように決めました。

n=1 のときを考えると m^2+2m< \mathrm{lim} より

2\leqq m<\sqrt{1+\mathrm{lim}}-1

n の範囲も m^2+2mn< \mathrm{lim} から決めました。

1\leqq n<\min\left\{n,\, \frac{\mathrm{lim}-m^2}{2m}\right\}

隣接リストを FixedArray で作る

隣接リストは mathematica version 12.1 で導入された FixedArray で作りました。

  adjacencyList = CreateDataStructure["FixedArray", {}, lim];

公式いわく「要素は素早く更新できる」らしいです。
たしかに adjacencyList = ConstantArray[ {}, lim] のように普通の多次元配列を使うのと比べるとちょっと速かったですね。

隣接リストへの要素の追加はまとめて行う

隣接リストに頂点を 1 つ 1 つ追加するかわりに (p,\, q) の倍数のリスト \left\{(p,\,q ),\, (2p,\,2q),\, \cdots\right\} をまとめて追加するようにしました。

  makeList[{m_, n_}] := Module[{p, q, lst},
    {p, q} = MinMax[{m^2 - n^2, 2 m*n + n^2}];
    lst = Times[{p, q}, #] & /@ Range@Quotient[lim, p + q];
    (adjacencyList["Part", #[[1]]] 
      = {adjacencyList["Part", #[[1]]], #[[2]]}) & /@ lst];

  makeList /@ mnTable;

最小の (p,\, q) は 3,\, 5) で，これから作られる (kp,\, kq) は 120{,}000/(3+5)=15{,}000 個あります。
これらの登録が 1 回で終わるようにしたわけです。

要素数 1 の隣接リストは調べない

p の隣接リストの要素数が 1 だったら「 p\to q, p\to r」をみたす q, r が存在しないので，こういう p は除外して考えました。

ちなみに要素数 1 の隣接リストは約 13,000 個で，要素数 2 個以上のものは約 22,000 個でした。

  vertex = Flatten@Position[adjacencyList, {Repeated[_, {2, Infinity}]}];

最終結果

以上をふまえて計算したところ約 2.5 秒でした。

  In[]:= ClearAll["Global`*"];
  AbsoluteTiming[
   lim = 12*10^4;
   mnTable = Flatten[Table[{m, n}, {m, 2, Sqrt[1 + lim] - 1},
      {n, 1, Min[m, Quotient[Quotient[lim, m] - m, 2]]}], 1];
   mnTable = 
    Select[mnTable, 
     CoprimeQ[#[[1]], #[[2]]] && ! Divisible[#[[1]] - #[[2]], 3] &];
   adjacencyList = CreateDataStructure["FixedArray", {}, lim];
   makeList[{m_, n_}] := Module[{p, q, lst},
     {p, q} = MinMax[{m^2 - n^2, 2 m*n + n^2}];
     lst = Times[{p, q}, #] & /@ Range@Quotient[lim, p + q];
     (adjacencyList[
          "Part", #[[1]]] = {adjacencyList["Part", #[[1]]], #[[2]]}) & /@
       lst];
   makeList /@ mnTable;
   adjacencyList = Union@Flatten[#, Infinity] & /@ Normal@adjacencyList;
   vertex = 
    Flatten@Position[adjacencyList, {Repeated[_, {2, Infinity}]}];
   calc[p_] := Module[{nums},
     nums = Select[Subsets[adjacencyList[[p]], {2}],
       MemberQ[adjacencyList[[#[[1]]]], #[[2]]] &];
     p + Total /@ nums];
   ans = Total[
     Select[DeleteDuplicates@Flatten[calc /@ vertex], # < lim &]];]
  
  Out[]= {2.51049, Null}

Project Euler 170 / 連結された積が0から9までのパンデジタル数になるパンデジタル数の最大値

mathematica project euler

問題概略

6 を 1273 と 9854 にかけると次のようになる。

\begin{align*}
&6 \times 1273 = 7638\\
&6 \times 9854 = 59124
\end{align*}

これらの積を連結すると，1 から 9 までを網羅するパンデジタル数 763859124 になる。
763859124 を「6 と (1273,\, 9854) の連結された積」と呼ぶことにする。


かける前の数を連結した 612739854 も 1 から 9 までを網羅するパンデジタル数である。
同じことが 0 から 9 までを網羅するパンデジタル数に対してもできる。


かける前の数を連結した数も 0 から 9 までのパンデジタル数となるような，ある整数と 2 つ以上の整数の連結された積が 0 から 9 までのパンデジタル数の最大値を求めよ。
https://projecteuler.net/problem=170

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

drive.google.com

右から攻めるか左から攻めるか

パンデジタル数を n 個に分割したときの各数を次のようにあらわします。

\begin{align*}
a\cdot b_1=c_1,\, a\cdot b_2=c_2,\, \cdots ,\, a\cdot b_n=c_n
\end{align*}

a と b_1～ b_n を連結した数を「左のパンデジタル数」とよんで， c_1～ c_n を連結した数を「右のパンデジタル数」とよぶことにします。
この問題で要求されているのは右のパンデジタル数の最大値です。

「左」→「右」と「右」→「左」のどちらで攻めるかを考えました。

「左」から「右」を求める方法：

  1. 左のパンデジタル数を a と b_1～ b_n に分割する
  2. a\cdot b_i\ (1\leqq i\leqq n) の結果を連結してパンデジタルかどうか調べる

「右」から「左」を求める方法：

  1. 右のパンデジタル数を c_1～ c_n に分割する
  2. c_1～ c_n の公約数 a を求める
  3. a と c_i/a=b_i\ (1\leqq i\leqq n) を連結してパンデジタルかどうか調べる
  4. 以上をすべての a に対して調べる

「左」→「右」の方が楽ですが，かけ算して連結しないと右のパンデジタル数がどんな数なのかわからないのが難点です。
本当に全探索しないと最大値は求められないので，計算量の面から「右」→「左」でいくことにしました。

実例を 1 つ探す

0～9 を使ったパンデジタル数は 9\cdot 9!=3,265,920 個あります。

これらに対してすべての分割法を試すのは多分無駄なので，「適当な仮定をおいて実例を 1 個探す」→「その数より大きい数についてちゃんと調べる」でいこうと思いました。

「分割数 n=2 かつ a は 2 桁以下」と仮定して計算しました。

  In[]:= AbsoluteTiming[
    lst = Select[Permutations@Reverse@Range[0, 9], First@# != 0 &];
    testk[lst_, k_] := Module[{c1, c2, aLst},
      {c1, c2} = FromDigits /@ TakeDrop[lst, k];
      aLst = Select[Divisors@GCD[c1, c2], # < 100 &];
      AnyTrue[aLst, 
       Sort[Join @@ IntegerDigits /@ {#, c1/#, c2/#}] == Range[0, 9] &]];
    isok[lst_] := AnyTrue[Range@9, testk[lst, #] &];
    ans = FromDigits@SelectFirst[lst, isok]]
   
  Out[]= {5.51237, 9857164023}

各行の内容はこうです。

  • 2 行目：解の候補となるパンデジタル数の桁数字のリストを作り，降順に並べる
  • 4 行目：桁数字のリストを先頭の k 個と残りにわけてそれぞれ数に直す（ c_1, c_2）
  • 5 行目： c_1, c_2 の公約数 a で 2 桁未満のもののリストを作る
  • 7 行目： a, c_1/a, c_2/a を連結してパンデジタルかどうか調べる
  • 9 行目：条件をみたす最初の桁数字のリストを抽出して数に直す

これでみつけた「9857164023」より大きなパンデジタル数は意外に多くて，10565 個ありました。
これらすべてに対して「右」→「左」のチェックをするのは大変なのでもう少しくわしく調べます。

分割数などをしぼりこむ

分割数をしぼりこむためのコードを書きました。
以下， x の桁数を len(x) であらわします。左のパンデジタル数の桁数の条件から

\begin{align*}
len(a)+\sum_{i=1}^n len(b_i)=10\mbox{　……(1)}]
\end{align*}

一般に x 桁の数と y 桁の数の積は x+y-1 桁か x+y 桁です（例： 2\cdot 1=2, 2\cdot 5=10）。

len(a\cdot b_i)=len(a)+len(b_i)-1 が k 回おこるとします。右のパンデジタル数の桁数の条件から

\begin{align*}
&\sum_{i=1}^n \left\{len(a)+len(b_i)\right\}-k=10\\
&\therefore n\cdot len(a)+\sum_{i=1}^n len(b_i)-k=10\mbox{　……(2)}
\end{align*}

(1)(2)を辺ごとに引くと (1-n)len(a)+k=0\mbox{　……(3)} を得ます。

2\leqq n\leqq 9, 1\leqq a\leqq 8, 0\leqq k\leqq n としてこれを解きました。

  eqn = {(1 - n) len_a + k == 0, 1 <= len_a <= 8, 
        2 <= n <= 9, 0 <= k <= n};
  sol = Reduce[eqn, {n, len_a, k}, Integers]

(3)をみたす \left(n,\, len(a),\, k\right) は次の 9 個でした。

\begin{align*}
&(2,\, 1,\, 1),\, (2,\, 2,\, 2),\, (3,\, 1,\, 2),\, (4,\, 1,\, 3),\, (5,\, 1,\, 4)\\
&(6,\, 1,\, 5),\, (7,\, 1,\, 6),\, (8,\, 1,\, 7),\, (9,\, 1,\, 8)
\end{align*}

(2,\, 1,\, 1) と 2,\, 2,\, 2) は「実例を 1 つ探す」で調べたケースです。

残り 7 パターンについては実は調べる必要がありません。

これらすべてが「 len(a)=1 かつ k=n-1」であることに注目しましょう。
n 個の a\cdot b_i=c_i のうち n-1 個では len(b_i)=len(c_i) が成り立ちます。

len(c_i)=len(a\cdot b_i)=len(a)+len(b_i)-1=len(b_i)

残り 1 個の a\cdot b_i=c_i では len(b_i)+1=len(c_i) が成り立ちます。

len(c_i)=len(a\cdot b_i)=len(a)+len(b_i)=len(b_i)+1

これらを使うと n\geqq 3 のケースはすべて n=2 のケースに帰着できることが示せます。

n=3 の場合を考えます。

a \cdot b_1 = c_1\mbox{　……(4)}
a \cdot b_2 = c_2\mbox{　……(5)}
a \cdot b_3 = c_3\mbox{　……(6)}

len(b_1)=len(c_1), len(b_2)=len(c_2), len(b_3)+1=len(c_3) として一般性を失いません。

len(b_2)=len(c_2)=L とおきます。(4)を 10^L 倍して(5)と辺ごとに足すと次のようになります。

a(b_1\cdot 10^L+b_2)=c_1\cdot 10^L+c_2

要するに 2 つの式を連結して1つにできるわけです。
これを繰り返すと n=4 ～ n=9 のケースはすべて「 n=2, len(a)=1, k=1」のケースに帰着できます。
結局，「実例を 1 つ探す」で求めた例がこの問題の解だったわけです。

高速化のための微調整

条件をみたすパンデジタル数を全部求めてみたら a はすべて 48 以下の 3 の倍数でした。

a=3,\, 6,\, 9,\, 12,\, 15,\, 18,\, 21,\, 24,\, 27,\, 36,\, 39,\, 42,\, 45,\, 48

これを証明して計算にいかします。

a は 3 の倍数

a\cdot b_1=c_1, a\cdot b_2=c_2 を辺ごとに足します。

a(b_1+b_2)=c_1+c_2\mbox{　……(7)}

パンデジタル数の条件から a, b_1, b_2 の桁数字の和も b_1, b_2 の桁数字の和も 0+1+\cdots+9=45 で 3 の倍数です。

a+b_1+b_2\equiv c_1+c_2\equiv 0\pmod{3}

(7)を \bmod\, 3 で考えます。

a(-a)\equiv 0\pmod{3}

a\equiv 0,\, \pm 1\pmod{3} のうちこれをみたすのは a\equiv 0\pmod{3} だけなので a は 3 の倍数です。

a は 50 未満

a\geqq 50 の場合を考えます。これは「分割数などをしぼりこむ」の \left(n,\, len(a),\, k\right)=(2,\, 2,\, 2) のケースです。

\begin{align*}
&a\cdot b_1=c_1,\, a\cdot b_2=c_2\\
&len(a\cdot b_1)=len(a)+len(b_1)-1=len(b_1)+1\mbox{　……(8)}\\
&len(a\cdot b_2)=len(a)+len(b_2)-1=len(b_2)+1
\end{align*}

b_1, b_2 の少なくとも 1 つは首位の数字が 2 以上です。
b_1 の首位の数字が 2 以上だとしましょう。
これを a\ (\geqq 50) 倍すると 2 桁以上増えます。

len(a\cdot b_1)\geqq len(b_1)+2

これは(8)と矛盾しているので a< 50 です。

再計算

以上をふまえてもう一度計算しました。計算時間は約 4 秒です。
もともとの計算時間が約 5.5 秒で，その半分は候補となるパンデジタル数の集合を作るのに費やされているので純粋な計算時間はだいたい半分になりました。

  In[]:= AbsoluteTiming[
    lst = Select[Permutations@Reverse@Range[0, 9], First@# != 0 &];
    testk[lst_, k_] := Module[{c1, c2, aLst},
      {c1, c2} = FromDigits /@ TakeDrop[lst, k];
      aLst = Intersection[Divisors@GCD[c1, c2], Range[3, 50, 3]];
      AnyTrue[aLst, 
       Sort[Join @@ IntegerDigits /@ {#, c1/#, c2/#}] == Range[0, 9] &]];
    isok[lst_] := AnyTrue[Range@9, testk[lst, #] &];
    ans = FromDigits@SelectFirst[lst, isok]]

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}