算額（その950）
一四 武州金鑽村 金鑽寺 文化10年(1813)
一五 武州金鑽村 金鑽寺 文化11年(1814)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
菱形の中に中央に4円（方併円），左右に2円（縦隔円），上下に2円（横隔円）を入れる。菱形の対角線の長い方が 60 寸，短い方が 32 寸のとき，それぞれの円の直径を求めよ。
菱形の対角線を 2a, 2b; a > b
方併円の半径と中心座標を r1, (r1, r1)
縦隔円の半径と中心座標を r2, (x2, 0)
横隔円の半径と中心座標を r3, (0, y3)
とおき，以下の連立方程式を解く。
eq1, (eq2, eq4), (eq3, eq5) は独立なので，別々に解くことができる。
include("julia-source.txt");
@syms a::positive, b::positive, r1::positive,
r2::positive, x2::positive,
r3::positive, y3::positive, c::positive
eq1 = dist2(a, 0, 0, b, r1, r1, r1)
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println
a/2 + b/2 - sqrt(a^2 + b^2)/2
# eq2 = dist2(a, 0, 0, b, x2, 0, r2)
eq2 = r2*sqrt(a^2 + b^2) - b*(a - x2)
eq4 = (x2 - r1)^2 + r1^2 - (r1 + r2)^2 |> expand
(ans_r2, ans_x2) = solve([eq2, eq4], (r2, x2))[2]
ans_r2 |> println
ans_x2 |> println
b*(a - (a^2*r1 - a*b^2 + b^2*r1 - b*r1*sqrt(a^2 + b^2) + sqrt(a^4*b^2 - 2*a^3*b^2*r1 + 2*a^3*b*r1*sqrt(a^2 + b^2) + a^2*b^4 + 2*a^2*b^2*r1^2 - 2*a^2*b*r1^2*sqrt(a^2 + b^2) - 2*a*b^4*r1 + 2*a*b^3*r1*sqrt(a^2 + b^2) + 2*b^4*r1^2 - 2*b^3*r1^2*sqrt(a^2 + b^2)))/a^2)/sqrt(a^2 + b^2)
(a^2*r1 - a*b^2 + b^2*r1 - b*r1*sqrt(a^2 + b^2) + sqrt(a^4*b^2 - 2*a^3*b^2*r1 + 2*a^3*b*r1*sqrt(a^2 + b^2) + a^2*b^4 + 2*a^2*b^2*r1^2 - 2*a^2*b*r1^2*sqrt(a^2 + b^2) - 2*a*b^4*r1 + 2*a*b^3*r1*sqrt(a^2 + b^2) + 2*b^4*r1^2 - 2*b^3*r1^2*sqrt(a^2 + b^2)))/a^2
# eq3 = dist2(a, 0, 0, b, 0, y3, r3)
eq3 = r3*sqrt(a^2 + b^2) - a*(b - y3)
eq5 = r1^2 + (y3 - r1)^2 - (r1 + r3)^2 |> expand;
(ans_r3, ans_y3) = solve([eq3, eq5], (r3, y3))[2]
ans_r3 |> println
ans_y3 |> println
a*(b - (-a^2*b + a^2*r1 - a*r1*sqrt(a^2 + b^2) + b^2*r1 + sqrt(a^4*b^2 - 2*a^4*b*r1 + 2*a^4*r1^2 + 2*a^3*b*r1*sqrt(a^2 + b^2) - 2*a^3*r1^2*sqrt(a^2 + b^2) + a^2*b^4 - 2*a^2*b^3*r1 + 2*a^2*b^2*r1^2 + 2*a*b^3*r1*sqrt(a^2 + b^2) - 2*a*b^2*r1^2*sqrt(a^2 + b^2)))/b^2)/sqrt(a^2 + b^2)
(-a^2*b + a^2*r1 - a*r1*sqrt(a^2 + b^2) + b^2*r1 + sqrt(a^4*b^2 - 2*a^4*b*r1 + 2*a^4*r1^2 + 2*a^3*b*r1*sqrt(a^2 + b^2) - 2*a^3*r1^2*sqrt(a^2 + b^2) + a^2*b^4 - 2*a^2*b^3*r1 + 2*a^2*b^2*r1^2 + 2*a*b^3*r1*sqrt(a^2 + b^2) - 2*a*b^2*r1^2*sqrt(a^2 + b^2)))/b^2
a, b が与えられたとき，r1, x2, r2, y3, r3 を順次求める式は以下のとおりである。
c = sqrt(a^2 + b^2)
r1 = a/2 + b/2 - c/2
x2 = (a^2*r1 - a*b^2 + b^2*r1 - b*r1*c + sqrt(a^4*b^2 - 2a^3*b^2*r1 + 2a^3*b*r1*c + a^2*b^4 + 2a^2*b^2*r1^2 - 2a^2*b*r1^2*c - 2a*b^4*r1 + 2a*b^3*r1*c + 2b^4*r1^2 - 2b^3*r1^2*c))/a^2
r2 = b*(a - x2)/c
y3 = (-a^2*b + a^2*r1 - a*r1*c + b^2*r1 + sqrt(a^4*b^2 - 2a^4*b*r1 + 2a^4*r1^2 + 2a^3*b*r1*c - 2a^3*r1^2*c + a^2*b^4 - 2a^2*b^3*r1 + 2a^2*b^2*r1^2 + 2a*b^3*r1*c - 2a*b^2*r1^2*c))/b^2
r3 = a*(b - y3)/c
a = 60/2, b = 32/2 のとき，中央の4個の円の直径（方併円径），左右にある2個の円の直径（縦隔円径），上下にある2個の円の直径（横隔円径）はそれぞれ，12 寸，12.5237 寸，5.91265 寸である。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (60, 32) ./2
c = sqrt(a^2 + b^2)
r1 = a/2 + b/2 - c/2
x2 = (a^2*r1 - a*b^2 + b^2*r1 - b*r1*c + sqrt(a^4*b^2 - 2a^3*b^2*r1 + 2a^3*b*r1*c + a^2*b^4 + 2a^2*b^2*r1^2 - 2a^2*b*r1^2*c - 2a*b^4*r1 + 2a*b^3*r1*c + 2b^4*r1^2 - 2b^3*r1^2*c))/a^2
r2 = b*(a - x2)/c
y3 = (-a^2*b + a^2*r1 - a*r1*c + b^2*r1 + sqrt(a^4*b^2 - 2a^4*b*r1 + 2a^4*r1^2 + 2a^3*b*r1*c - 2a^3*r1^2*c + a^2*b^4 - 2a^2*b^3*r1 + 2a^2*b^2*r1^2 + 2a*b^3*r1*c - 2a*b^2*r1^2*c))/b^2
r3 = a*(b - y3)/c
@printf("方併円径 = %g; 縦隔円径 = %g; 横隔円径 = %g\n", 2r1, 2r2, 2r3)
@printf("r1 = %g; r2 = %g; x2 = %g; r3 = %g; y3 = %g\n", r1, r2, x2, r3, y3)
plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:black, lw=0.5)
circle4(r1, r1, r1)
circle2(x2, 0, r2, :green)
circle22(0, y3, r3, :blue)
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(a, 0, " a", :black, :left, :bottom, delta=delta/2)
point(0, b, " b", :black, :left, :bottom, delta=delta/2)
point(r1, r1, "方併円:r1\n(r1,r1)", :red, :center, delta=-delta)
point(x2, 0, "縦隔円:r2\n(x2,0)", :green, :center, delta=-delta)
point(0, y3, " 横隔円:r3,(0,y3)", :blue, :left, :vcenter)
end
end;
