算額（その812）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
長方形内に大円，中円，小円が入っている。長方形の長辺，短辺がそれぞれ 86634 寸，77008 寸のとき，大円の直径はいかほどか。
長方形の長辺，短辺を a, b
大円の半径と中心座標を r1, (a - r1, r1)
中円の半径と中心座標を r2, (r2, b - r2)
小円の半径と中心座標を r3, (r3, r3)
とおいて以下の連立方程式を解く。
なお，SymPy の能力的に，a, b を変数のままにして解くことはできない。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, r2::positive, r3::positive,
a::positive, b::positive
(a, b) = (86634, 77008)
eq1 = (a - r1 - r2)^2 + (b - r2 - r1)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (r2 - r3)^2 + (b - r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))
res[2]
(-57756*sqrt(338 + 240*sqrt(2)) + 24065 + 81821*sqrt(169 + 120*sqrt(2)), -81821*sqrt(169 + 120*sqrt(2)) + 24065 + 57756*sqrt(338 + 240*sqrt(2)), -4813*sqrt(169 + 120*sqrt(2)) + 24065 + 57756*sqrt(2))
6 通りの解が得られるが，2 番目のものが適解である。
長方形の長辺，短辺がそれぞれ 86634 寸，77008 寸のとき，大円の直径は 53345.0001130072 寸である。
2res[2][1].evalf() |> println
53345.0001130072
その他のパラメータは以下のとおりである。
r1 = 26672.5; r2 = 21457.5; r3 = 17166.1
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (86634, 77008)
t = sqrt(169 + 120√2)
u = 81821 - 57756√2
(r1, r2, r3) = 24065 .+ (t*u, -t*u, 57756√2 - 4813t)
@printf("大円の直径 = %g; 中円の直径 = %g; 小円の直径 = %g\n", 2r1, 2r2, 2r3)
@printf("r1 = %g; r2 = %g; r3 = %g\n", r1, r2, r3)
plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
circle(a - r1, r1, r1)
circle(r2, b - r2, r2, :blue)
circle(r3, r3, r3, :green)
if more == true
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(a - r1, r1, "大円:r1,(a-r1,r1)", :red, :center, delta=-delta)
point(r1, b - r2, "中円:r2,(r2,b-r2)", :blue, :center, delta=-delta)
point(r3, r3, "小円:r3,(r3,r3)", :green, :center, delta=-delta)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
point(0, b, " b", :black, :left, :bottom, delta=delta)
end
end;
