算額（その1364）
三十四 群馬県多野郡新町 稲荷神社 文政3年(1820)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円15個，等脚台形
#Julia, #SymPy, #算額, #和算
等脚台形の中に，大円 5 個，小円 3 個，甲円 2 個，乙円 5 個を容れる。台形の上底が 14 寸，下底が 32 寸のとき，小円の直径はいかほどか。
注：「問」では台形の上底と下底が与えられているが，右側の 3 個の大円の中心を結ぶと正三角形になり，台形の斜辺と下底がなす角は 60° なので，下底が 32 寸のとき，上底は 14 寸にはならない。そこで，下底のみが与えられるとして，問を解く。算額（その1366）も参照のこと。
等脚台形の下底，上底，高さを 2a, 2b, h
大円の半径と中心座標を r1, (2r1, r1), (r1, h - r1)
小円の半径と中心座標を r2, (0, 2r1 + r2)
甲円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (0, h - r4)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, h::positive,
r1::positive, r2::positive,
r3::positive, x3::positive, r4::positive;
eq1 = r1^2 + (h - 2r1)^2 - 4r1^2
eq2 = r1/(a - 2r1) - 1/√Sym(3)
eq3 = r1/(b - r1) - √Sym(3);
# res1 = solve([eq1, eq2, eq3], (r1, b, h));
eq4 = r1^2 + (h - r1 - (2r1 + r2))^2 - (r1 + r2)^2
eq5 = r1^2 + (r1 - r4)^2 - (r1 + r4)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r4, b, h))
eq6 = r3/(a - x3) - 1/√Sym(3)
eq7 = x3 - 2r1 - 2sqrt(r1*r3)
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7],
(r1, r2, r3, x3, r4, b, h))[2] # 2 of 2
(a*(2 - sqrt(3)), -1129438*a*sqrt(79315912984*sqrt(3) + 137379191137) - 2*a + sqrt(3)*a + 1956244*sqrt(3)*a*sqrt(79315912984*sqrt(3) + 137379191137)/3, a*(-5*sqrt(3) - 4*sqrt(7 - 4*sqrt(3)) + 10)/3, 2*a*(-5*sqrt(3) + 2*sqrt(21 - 12*sqrt(3)) + 9)/3, a*(2 - sqrt(3))/4, a*(3 - sqrt(3))/3, a*(-1694157*sqrt(79315912984*sqrt(3) + 137379191137) - 2*sqrt(3) + 4 + 978122*sqrt(412137573411 + 237947738952*sqrt(3))))
パラメータによっては二重根号を外して簡約化できる。
r1
res[1] |> println
a*(2 - sqrt(3))
r2
res[2] |> sympy.sqrtdenest |> simplify |> println
a*(-12 + 7*sqrt(3))/3
r3
res[3] |> sympy.sqrtdenest |> simplify |> println
a*(2 - sqrt(3))/3
x3
res[4] |> sympy.sqrtdenest |> simplify |> println
2*a*(3 - sqrt(3))/3
r4
res[5] |> println
a*(2 - sqrt(3))/4
b
res[6] |> println
a*(3 - sqrt(3))/3
h
res[7] |> sympy.sqrtdenest |> simplify |> println
a
小円の半径は a*(7√3 - 12)/3 なので，下底が 32 寸のとき小円の直径は 32*(7√3 - 12)/3 = 1.3264602984761684 寸である。
このとき，上底は 32*(3 - sqrt(3))/3 = 13.524791385931977 である。
全てのパラメータは以下のとおりである。
a = 16; r1 = 4.28719; r2 = 0.66323; r3 = 1.42906; x3 = 13.5248; r4 = 1.0718; b = 6.7624; h = 16
---
なお，小円の直径を知るだけならば，甲円，乙円は不要である（問題を複雑そうに見せる虚仮威しである）。
以下のような手順で小円の直径を（紙と鉛筆だけで）求めることができる。
このことからも，問で与えられた上底の条件は無用なものどころか，有害なものであることがわかる。
下底 2a と大円の半径 r1 の関係式は以下であり，大円の半径は r1 = a*(2 - √3) である。
@syms a, r1
eq01 = (2 + √Sym(3))r1 - a
eq01 |> println
-a + r1*(sqrt(3) + 2)
r1 = solve(eq01, r1)[1] |> simplify
r1 |> println
a*(2 - sqrt(3))
一辺の長さが 2r1 の正三角形の頂点を中心とする3個の大円が作る中央の隙間にピッタリハマる円の半径は以下のようになる（図中の灰色で描いた小さな直角三角形を参照）。
r1*2/√Sym(3) - r1 |> simplify |> println
a*(-12 + 7*sqrt(3))/3
下底が 32 寸のとき，小円の直径は 1.3264602984761684 寸である。
a = 32/2
2*a*(-12 + 7*sqrt(3))/3
1.3264602984761684
以下は図を描くプログラム
function draw(a, more)
pyplot(size=(800, 800), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3, x3, r4, b, h) = (a*(2 - sqrt(3)), -1129438*a*sqrt(79315912984*sqrt(3) + 137379191137) - 2*a + sqrt(3)*a + 1956244*sqrt(3)*a*sqrt(79315912984*sqrt(3) + 137379191137)/3, a*(-5*sqrt(3) - 4*sqrt(7 - 4*sqrt(3)) + 10)/3, 2*a*(-5*sqrt(3) + 2*sqrt(21 - 12*sqrt(3)) + 9)/3, a*(2 - sqrt(3))/4, a*(3 - sqrt(3))/3, a*(-1694157*sqrt(79315912984*sqrt(3) + 137379191137) - 2*sqrt(3) + 4 + 978122*sqrt(412137573411 + 237947738952*sqrt(3))))
r1 = a*(2 - √3)
r2 = a*(7√3 - 12)/3
r3 = a*(2 - √3)/3
x3 = 2*a*(3 - √3)/3
r4 = a*(2 - √3)/4
b = a*(3 - √3)/3
h = a
@printf("下底が %g のとき，小円の直径は %g である（上底は %g）。\n", 2a, 2r2, 2b)
@printf("a = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g; r4 = %g; b = %g; h = %g\n", a, r1, r2, r3, x3, r4, b, h)
plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:brown, lw=0.5)
circle(0, r1, r1)
circle2(2r1, r1, r1)
circle2(r1, h - r1, r1)
circle(0, 2r1 + r2, r2, :blue)
circle2(r1, h - 2r1 - r2, r2, :blue)
circle(0, h - r4, r4, :green)
circle2(r1, r4, r4, :green)
circle2(x3, r3, r3, :orange)
l = h - (2r1 + r2) - r4
circle2(r1 + l*cosd(30), h - (2r1 + r2) + l*sind(30), r4, :green)
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(2r1, r1, "大円:r1,(2r1,r1)", :red, :center, delta=-delta)
point(0, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta)
point(r1, (1 + √3)r1, "大円:r1,(r1,(1+√3)r1)", :red, :center, delta=-delta)
point(0, 2r1 + r2, "小円:r2,(0,2r1+r2)", :blue, :left, :vcenter, deltax=4delta)
point(0, h - r4, "乙円:(0,h-r4)", :green, :left, :vcenter, deltax=7delta)
point(x3, r3, "甲円:(x3,r3)", :orange, :right, :vcenter, deltax=-9delta)
point(b, h, "(b,h)", :brown, :right, :bottom, delta=delta/2)
point(a, 0, "a", :brown, :left, :bottom, delta=delta/2)
plot!([r1, 2r1, r1, r1], [r1, r1, h - 2r1 - r2, r1], color=:gray70, lw=0.5)
end
end;
draw(32/2, true)
