算額（その858）
二十二 岩手県一関市瑞山 駒形根神社 明治41年(1908)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
正三角形内に甲円，乙円，丙円を各 1 個，丁円を 2 個入れる。甲円の直径がわかっているときに丁円の直径を求めよ。
注：引用元の図は誤解を招く不正確なものである。以下のように丙円，丁円は外接している。「術」に述べられている解になるためには，このような図でなければならない。
正三角形の一辺の長さを 2a
甲円の半径と中心座標を r1, (0, 2r3 + 2r2 + r1)
乙円の半径と中心座標を r2, (0, 2r3 + r2)
丙円の半径と中心座標を r3, (0, r3)
丁円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive,
r4::positive, x4::positive, d
s3 = √Sym(3)
eq1 = 2r1 + r1 + 2r2 + 2r3 - s3*a
eq2 = (s3*a - (2r3 + r2)) - 2r2
eq3 = dist(a, 0, 0, s3*a, x4, r4) - r4^2 |> (x -> apart(x, d)) |> numerator
eq4 = x4^2 + ((2r3 + r2) - r4)^2 - (r2 + r4)^2
eq5 = x4^2 + (r4 - r3)^2 - (r3 + r4)^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r2, r3, r4, x4));
2 組の解が得られるが，2 番目のものが適解である。そのうちの 4 番目が丁円の半径。
res[2] |> println
res[2][4] |> println
(r1*(-6*sqrt(265 + 156*sqrt(3)) + sqrt(795 + 468*sqrt(3)) + 132 + 110*sqrt(3))/33, 3*r1, r1*(-2*sqrt(795 + 468*sqrt(3)) + 11 + sqrt(265 + 156*sqrt(3)) + 44*sqrt(3))/22, r1*(-10*sqrt(265 + 156*sqrt(3)) - 2*sqrt(795 + 468*sqrt(3)) + 187 + 110*sqrt(3))/33, -3*sqrt(3)*r1*(-20*sqrt(265/324 + 13*sqrt(3)/27) + 70/9 + 20*sqrt(3)/3)/10)
r1*(-10*sqrt(265 + 156*sqrt(3)) - 2*sqrt(795 + 468*sqrt(3)) + 187 + 110*sqrt(3))/33
甲円の半径に対する丁円の半径の倍数は，以下のように簡約化できる。
res1 = res[2][4]/r1
res1 |> println
-10*sqrt(265 + 156*sqrt(3))/33 - 2*sqrt(795 + 468*sqrt(3))/33 + 17/3 + 10*sqrt(3)/3
res2 = (res1.args[2] + res1.args[3])^2 |> simplify |> sqrt |> simplify
res2 |> println
2*sqrt(100 + 58*sqrt(3))/3
res3 = res1.args[1] + res1.args[4] - res2
res3 |> println
-2*sqrt(100 + 58*sqrt(3))/3 + 17/3 + 10*sqrt(3)/3
甲円の半径に対する丁円の半径の倍数は約 2 である。正確に 2 ではないことに驚く。まさに奇跡。
res3 |> N
2.001267660707633
術には以下の式が書かれている。見かけは簡単であるが，上述の結果と同じである。
天 = √3 + 1
(sqrt(8天 + 5) - 天)^2/3
2.001267660707633
function draw(more=false)
pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 0.5
s3 = √3
s = sqrt(795 + 468s3)
t = sqrt(265 + 156s3)
(a, r2, r3, r4, x4) = r1 .* (
(s - 6t + 132 + 110s3)/33,
3,
(t - 2s + 11 +44s3)/22,
(187 + 110s3 - 10t - 2s)/33,
-3s3*(70/9 + 20s3/3 - 20sqrt(265/324 + 13s3/27))/10)
@printf("甲円の直径が %g のとき，丁円の直径は %g\n", 2r1, 2r4)
@printf("r1 = %g; a = %g; r2 = %g; r3 = %g; r4 = %g; x4 = %g\n", r1, a, r2, r3, r4, x4)
plot([a, 0, -a, a], [0, s3*a, 0, 0], color=:blue, lw=0.5)
circle(0, s3*a - 2r1, r1)
circle(0, s3*a - 3r1 - r2, r2, :green)
circle(0, r3, r3, :orange)
circle2(x4, r4, r4, :magenta)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(0, s3*a - 2r1, " 甲円:r1,(0,√3*a-2r1)", :red, :left, :vcenter)
point(0, s3*a - 3r1 - r2, "乙円:r2,(0,√3*a-3r1-r2)", :green, :center, delta=-delta/2)
point(0, r3, "丙円;r3\n(0,r3)", :orange, :center, delta=-delta/2)
point(x4, r4, "丁円;r4,(x4,r4)", :magenta, :center, delta=-delta/2)
point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
point(0, √3a, " √3a", :blue, :left, :vcenter)
end
end;
