算額（その614）

算額（その614）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/7

2 本の直線に挟まれて大円，中円，小円が互いに接している。中円と小円の直径が与えられたとき，大円の直径を求めよ。

時計回りに 90° 回転させた図で考える。
大円の半径と中心座標を r1, (x1, 0)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (x3, 0)
2 直線の交点座標を (a, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, x1::positive, r2::positive,
     r3::positive, x3::positive
(r2, r3) = (5, 2)
tt = r3/(a - x3)
eq1 = r2/a - r3/(a - x3)
eq2 = x1^2 + r2^2 - (r1 + r2)^2
eq3 = (x3 - x1)^2 + r3^2 - (r1 + r3)^2
eq4 = r1/sqrt((a - x1)^2 - r1^2) - 2tt/(1 - tt^2)
res = solve([eq1, eq2, eq3, eq4], (a, r1, x1, x3))

   1-element Vector{NTuple{4, Sym}}:
    (5*sqrt(169 + 56*sqrt(10))/3, -21/8 - 49*sqrt(10)/160 + 7*sqrt(10)*sqrt(56*sqrt(10) + 209)/160 + 3*sqrt(56*sqrt(10) + 209)/8, 3*sqrt(169 + 56*sqrt(10))*sqrt(56*sqrt(10) + 209)/(2*(160 + 56*sqrt(10))) - (139 + 56*sqrt(10))*sqrt(169 + 56*sqrt(10))/(2*(3 - sqrt(169 + 56*sqrt(10)))*(3 + sqrt(169 + 56*sqrt(10)))), sqrt(169 + 56*sqrt(10)))

大円の半径の式

res[1][2] |> println

   -21/8 - 49*sqrt(10)/160 + 7*sqrt(10)*sqrt(56*sqrt(10) + 209)/160 + 3*sqrt(56*sqrt(10) + 209)/8

簡約化できる

res[1][2] |> sympy.sqrtdenest |> simplify |> println

   7/4 + 3*sqrt(10)/2

r2,r3 を未知数として解くと不適切解しか得られない。
r2, r3 として整数を与えて解を吟味することにより一般解を得ることができる。
r2 = 9, r3 = 1〜7 を試せば，r1 は (r2 + r3)/4 + 3sqrt(r2*r3)/2 らしいことがわかる。
また，整数値でなくても成り立つこともわかる。

術では「中径と小径を乗じて得られる数を開平して二倍し，これを極と名付ける。この値から中径，小径を引き四で割る。極からこの値を引けば大円の径になる」と記している。
Julia で書くと
極 = 2sqrt(r2*r3)
極 - (極 - r3 - r2)/4
のようになる。これは十分簡約化された数式であるが，Julia はもっと簡約化する。すなわち r2/4 + r3/4 + 3*sqrt(r2*r3)/2 でこれは前述した式に等しい。

@syms r2, r3
極 = 2sqrt(r2*r3)
極 - (極 - r3 - r2)/4 |> println

   r2/4 + r3/4 + 3*sqrt(r2*r3)/2

(r2/4 + r3/4 + 3*sqrt(r2*r3)/2)(r2 =>5, r3=>2) |> println

   7/4 + 3*sqrt(10)/2

中円，小円の直径が 10, 4 のとき，大円の直径は 12.9868。

その他のパラメータは以下の通り。

   a = 31.0057;  r1 = 6.49342;  x1 = 10.3488;  x3 = 18.6034

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (5, 2)
   (a, r1, x1, x3) = res[1]
   @printf("a = %g;  r1 = %g;  x1 = %g;  x3 = %g\n", a, r1, x1, x3)
   plot()
   circle(0, r2, r2)
   circle(0, -r2, r2)
   circle(x1, 0, r1, :blue)
   circle(x3, r3, r3, :green)
   circle(x3, -r3, r3, :green)
   abline(a, 0, r1/sqrt((a - x1)^2 - r1^2), -r2, a)
   abline(a, 0, -r1/sqrt((a - x1)^2 - r1^2), -r2, a)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, 0, " 大円:r1,(r1,0)", :blue, :center, :bottom, delta=delta)
       point(0, r2, " 中円:r2\n (0,r2)", :red, :left, :vcenter)
       point(x3, r3, "小円:r3,(x3,r3)", :black, :left, delta=-2delta)
       point(a, 0, "a", :black, :left, :bottom, delta=delta)
   end
end;