算額（その1145）
一〇一 大宮市高鼻町 氷川神社 明治31年(1898)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円2個，楕円2個，長方形，斜線2本
長方形の中に2本の斜線を引き，大円 2 個と楕円 2 個を容れる。楕円の長径が 4 寸，大円の直径が 8 寸のとき，楕円の短径はいかほどか。
長方形の長辺の長さを 2c, （短辺の長さは 2r）
斜線と長方形の長辺との交点の座標を (d, 0)
大円の半径と中心座標を r, (c - r, r)
楕円の長半径，短半径，中心座標を a, b, (0, b), (0, 2b + a)
楕円と斜線の接点座標を (x1, y1), (x2, y2)
とおき，以下の連立方程式を解く。
include("julia-source.txt")
using SymPy
@syms a::positive, b::positive, r::positivr,
c::positive, d::positive,
x1::positive, y1::positive,
x2::positive, y2::positive
@syms a, b, r, c, d, x1, y1, x2, y2
slope = -2r/d
eq1 = x1^2/a^2 + (y1 - b)^2/b^2 - 1
eq2 = -b^2*x1/(a^2*(y1 - b)) - slope
eq3 = y1/(d - x1) + slope
eq4 = x2^2/b^2 + (y2 - a - 2b)^2/a^2 - 1
eq5 = -a^2*x2/(b^2*(y2 - a - 2b)) - slope
eq6 = y2/(d - x2) + slope
eq7 = dist2(0, 2r, d, 0, c - r, r, r);
まず，eq1, eq2, eq3 を解き，b, x1, x2 を求める。
res1 = solve([eq1, eq2, eq3], (b, x1, y1))[1]
(-a^2*r/d^2 + r, 2*a^2*d/(a^2 + d^2), 2*r*(-a^2 + d^2)/(a^2 + d^2))
eq4, eq5 に b を代入し，eq6, eq7 とともに連立方程式を解き，x2, y2, c, d を求める。
eq4 = eq4(b => res1[1]);
eq5 = eq5(b => res1[1]);
res2 = solve([eq4, eq5, eq6, eq7], (x2, y2, c, d))[8] # 8 of 8
(2*sqrt(a^2*r^2/(a^2 - a*r + r^2))*(a^4 - 3*a^3*r + 4*a^2*r^2 - 3*a*r^3 + r^4)/(r^2*(2*a^2 - 3*a*r + 2*r^2)), 2*a*(-2*a^3 + 6*a^2*r - 6*a*r^2 + 3*r^3)/(r*(2*a^2 - 3*a*r + 2*r^2)), r*sqrt((5*a^2 - 4*a*r + 4*r^2)/(a^2 - a*r + r^2))/2 + r + sqrt(a^2*r^2/(a^2 - a*r + r^2))/2, sqrt(a^2*r^2/(a^2 - a*r + r^2)))
短半径は最初の連立方程式の解で -a^2*r/d^2 + r となっており，まだ未知数 d を含んでいる。
d は次の連立方程式の解で sqrt(a^2*r^2/(a^2 - a*r + r^2)) であることがわかる。
そこで，d を短半径の式に代入すれば a, r のみを含む式 a*(r - a)/r = (a*r - a^2)/r になる。これは「術」と同じである。
長半径が a = 4/2，大円の半径が r = 8/2 のとき，短半径は a*(-a + r)/r = 1 である（短径は 2）。
res1[1](d => res2[4]) |> simplify |> println
a*(-a + r)/r
術は以下のようになっている。
長径 = 4
大径 = 8
(長径*大径 - 長径^2)/大径
2.0
その他のパラメータは以下の通りである。
a = 2; r = 4; b = 1; x1 = 1.97949; y1 = 1.14286; x2 = 0.866025; y2 = 5; c = 9.31803; d = 2.3094
function draw(a, r, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x2, y2, c, d) = (2*sqrt(a^2*r^2/(a^2 - a*r + r^2))*(a^4 - 3*a^3*r + 4*a^2*r^2 - 3*a*r^3 + r^4)/(r^2*(2*a^2 - 3*a*r + 2*r^2)), 2*a*(-2*a^3 + 6*a^2*r - 6*a*r^2 + 3*r^3)/(r*(2*a^2 - 3*a*r + 2*r^2)), r*sqrt((5*a^2 - 4*a*r + 4*r^2)/(a^2 - a*r + r^2))/2 + r + sqrt(a^2*r^2/(a^2 - a*r + r^2))/2, sqrt(a^2*r^2/(a^2 - a*r + r^2)))
(b, x1, y1) = (-a^2*r/d^2 + r, 2*a^2*d/(a^2 + d^2), 2*r*(-a^2 + d^2)/(a^2 + d^2))
@printf("楕円の長径が %g，大円の直径が %g のとき，楕円の短径は %g である。\n", 2a, 2r, 2b)
@printf("a = %g; r = %g; b = %g; x1 = %g; y1 = %g; x2 = %g; y2 = %g; c = %g; d = %g\n", a, r, b, x1, y1, x2, y2, c, d)
plot([c, c, -c, -c, c], [0, 2r, 2r, 0, 0], color=:blue, lw=0.5)
plot!([-d, 0, d], [0, 2r, 0], color=:magenta, lw=0.5)
circle2(c - r, r, r)
ellipse(0, b, a, b, color=:green)
ellipse(0, a + 2b, b, a, color=: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(c, 0, " c", :blue, :left, :bottom, delta=delta/2)
point(d, 0, " d", :magenta, :left, :bottom, delta=delta/2)
point(c - r, r, "大円:r,(c-r,r)", :red, :center, delta=-delta/2)
point(0, 2b, "2b", :green, :center, delta=-2delta)
point(0, 2b + 2a, "2b+2a", :black, :center, :bottom, delta=delta)
end
end;
