Project Euler 199 / 反復円充填

問題概略
パターンを見抜く
デカルトの円定理
計算に使う関係式
python と mathematicaで実装
- python で解く
- mathematicaで解く

問題概略

半径が等しくて互いに外接する 3 つの円が大きい円に内接している。
この図形がもつ 4 つのすき間を円で埋める操作を繰り返す。
各ステップですべてのすき間にそれぞれ最大の円を置いていくと，3 ステップ後には図のようになる。
108 のすき間ができ，円で埋められていない面積の比率は 10 進 8 桁に四捨五入して 0.06790342 である。
10 ステップ後に円で埋まっていない面積の比率を 10 進 8 桁に四捨五入して $x.xxxxxxxx$ の形で答えよ。
https://projecteuler.net/problem=199

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

drive.google.com

パターンを見抜く

この問題で求めるのは比率だけなので，一番外側にある大円の半径を $r_0=1$ として一般性を失いません。
内部にある 3 つの小円の半径は $r_1=2\sqrt{3}-3$ です。

与えられた図を漫然と眺めるととりつく島もない感じですが，色ごとに見るとパターンがわかります。
「3 つの円から新しい円ができる」→「その円のまわりに新しい円が 3 つできる」の繰り返しです。

最初にある 4 つの円から作られる円の半径をまとめて $r$ であらわすと，求める比率は次の式で与えられます。

\begin{align*}
\texttt{ans}=\frac{\pi{r_0}^2-3\pi{r_1}^2-\sum \pi r^2 }{\pi{r_0}^2} &=1-3{r_1}^2-\sum r^2
\end{align*}

デカルトの円定理

新しくできる円の半径はデカルトの円定理で求められます。

w.wiki

これは「互いに接する 4 つの円の半径はある 2 次方程式をみたす」という定理です。

半径 $r$ の円の曲率を $k=\pm \frac{1}{r}$ で定義します。
符号は他の円を内部に含むときマイナスで，含まないときはプラスです。
この問題の場合，一番外側にある大円だけマイナスで，残りの円はすべてプラスです。

$k_0=-\frac{1}{r_0}=-1,\, k_i=\frac{1}{r_i}\quad (i\ne 0)$

曲率 $a$ , $b$ , $c$ , $k$ の円が互いに接するとき

$(a+b+c+k)^2=2(a^2+b^2+c^2+k^2)$

が成立します。正の解を求めると

$k=a+b+c+2\sqrt{ab+bc+ca}$

新しくできる円の半径は $r=1/k$ です。

計算に使う関係式

$k$ と $a$ , $b$ , $c$ の関係を $k=f(a,\, b,\, c)$ であらわすと，次のようにして円が作られていきます。

曲率 $a$ , $b$ , $c$ の円から曲率 $k=f(a,\, b,\, c)$ の円ができる
曲率 $f(k,\, a,\, b)$ , $f(k,\, b,\, c)$ , $f(k,\, c,\, a)$ の円ができる

求める比はこれを 10 回繰り返したときの $1-3{r_1}^2-\sum \frac{1}{k^2}$ です。
dfs と対称性を使って計算しました。

上の操作を $a$ , $b$ , $c$ からはじめて $n$ 回繰り返したときの $r^2=1/k^2$ の和を $f(a,\, b,\, c,\, n)$ で表します。

$\texttt{ans}=1-3{r_1}^2-3f(k_0,\, k_1,\, k_1,\, 10)-f(k_1,\, k_1,\, k_1,\, 10)$

python と mathematicaで実装

python で解く

pythonによるコードがこちら。計算時間は 0.03 秒～0.04 秒でした。

    import time, math

    def f(a, b, c, depth):
        if depth == 0:
            return 0
        k = a + b + c + 2 * math.sqrt(a * b + b * c + c * a)
        return 1 / (k * k) + f(k, a, b, depth - 1) \
               + f(k, b, c, depth - 1) + f(k, c, a, depth - 1)
    
    def solve(n):
        r0, r1 = 1, 2 * math.sqrt(3) - 3
        k0, k1 = -1 / r0, 1 / r1
        return 1 - 3 * r1 ** 2 - 3 * f(k0, k1, k1, n) - f(k1, k1, k1, n)
    
    if __name__ == '__main__':
        start = time.time()
        print('{:.8f}'.format(solve(10)))
        elapsed_time = time.time() - start
        print("elapsed_time: {0} [sec]".format(elapsed_time))

mathematicaで解く

mathematica は厳密解を求めようとするようで，何も考えずに計算させるとすごい時間がかかります。
円を作る操作が 8 回を超えたあたりで固まったように遅くなりました。

$r_1$ の値を N で囲んで，数値計算するように指定したら 0.5 秒くらいで答えが出ました。

    r0 = 1;
    r1 = N[2 Sqrt@3 - 3, 9];
    k0 = -1/r0;
    k1 = 1/r1;
    f[a_, b_, c_, 0] := 0;
    f[a_, b_, c_, n_] := Module[{k = a + b + c + 2*Sqrt[a*b + b*c + c*a]},
       1/k^2 + f[k, a, b, n - 1] + f[k, b, c, n - 1] + f[k, c, a, n - 1]];
    n = 10;
    ans = N[1 - 3*r1^2 - 3*f[k0, k1, k1, n] - f[k1, k1, k1, n], 8]

砂場で遊ぼう

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